############################################################################
# Bypassing the Enemy (2017)
# This file contains the replication code for the paper 
# "Bypassing the Enemy: Distributive Politics, Credit Claiming, and Nonstate Organizations in Brazil"

############ Getting data

########### Recoding variables

########### Tables and Figures

########### Data and analyses cited in the appendix and main paper

########### Additional analyses not shown in appendix

############################################################################

rm(list = ls(all = TRUE))

#Preambule
#R version 3.3.3

#Change to the appropiate directory
dir <- ""

options(scipen=999) # supressing scientific notation
par(mar=c(5.1,4.1,4.1,2.1)) 
par(mfrow=c(1,1))

#Updates to version 2

#Use the appropriate package versions (by either installing them manually or checkpoint package)
#Producing pdf figures instead of publication eps because of package conflicts (when using checkpoint)
#If wish to produce publication-quality figures, use the Cairo package and cairo_ps instead of pdf


#If you're getting the warning below, it's because you're using
#a R version other than 3.3.2
#Warning messages:
#  1: In (X - c)/h :
#  Recycling array of length 1 in vector-array arithmetic is deprecated.
#Use c() or as.vector() instead.

#You may need to setwd(dir) to get checkpoint to work properly

#setwd(dir)
library(checkpoint)
checkpoint("2017-03-30", checkpointLocation = tempdir())

#Run these packages
library(sandwich) 
library(xtable)
library(ri)
library(gtools)
library(rdrobust)
library(Hmisc)
library(rdd)
library(texreg)
library(memisc)
library(rgenoud)
library(stargazer)
library(ebal)
library(cjoint)
library(gridExtra)
library(tidyverse)
library(plotrix)

#Functions
source(paste0(dir, "2.bypassing_functions_final.R"))

############################################################################
#Getting main data
load(paste0(dir, "final_replication.RData"))

#Getting auxiliary data
load(paste0(dir, "final_replication_auxiliary.RData"))

############################################################################
#Creating variables used in the analysis (federal and state transfers)

#Windows (used in plots and tables)
H <- c(seq(from = .005, to = .0095, by = .0005), seq(from = .01, to = .05, by = .001))

#Clusters
data_200315$ibge2 <- as.numeric(data_200315$ibge)
data_200715$ibge2 <- as.numeric(data_200715$ibge)

#Nonzero for outcome variable
data_200315$nonzero <- ifelse(data_200315$dtransfers.pc==0, 0, 1)
data_200315$nonzero.mayors <- ifelse(data_200315$dtransfers.pc.mayors==0, 0, 1)
data_200715$nonzero <- ifelse(data_200715$dtransfers.pc==0, 0, 1)
data_200715$nonzero.mayors <- ifelse(data_200715$dtransfers.pc.mayors==0, 0, 1)

#Recoding running variable
data_200315$vote_margin_share2 <- data_200315$vote_margin_share*100
data_200715$vote_margin_share2 <- data_200715$vote_margin_share*100

#Recoding for per capita or vote share

data_200315$p_13.1 <- data_200315$p_13/data_200315$p_val
data_200315$p_45.1 <- data_200315$p_45/data_200315$p_val
data_200315$f.nom_13.1 <- data_200315$f.nom_13/data_200315$f.nom_total
data_200315$f.nom_45.1 <- data_200315$f.nom_45/data_200315$f.nom_total
data_200315$e.nom_13.1 <- data_200315$e.nom_13/data_200315$e.nom_total
data_200315$e.nom_45.1 <- data_200315$e.nom_45/data_200315$e.nom_total
data_200315$g_13.1 <- data_200315$g_13/data_200315$g_val
data_200315$g_45.1 <- data_200315$g_45/data_200315$g_val
data_200315$Comparecimento1t.1 <- data_200315$Comparecimento1t/data_200315$Aptos

data_200715$p_13.1 <- data_200715$p_13/data_200715$p_val
data_200715$p_45.1 <- data_200715$p_45/data_200715$p_val
data_200715$f.nom_13.1 <- data_200715$f.nom_13/data_200715$f.nom_total
data_200715$f.nom_45.1 <- data_200715$f.nom_45/data_200715$f.nom_total
data_200715$e.nom_13.1 <- data_200715$e.nom_13/data_200715$e.nom_total
data_200715$e.nom_45.1 <- data_200715$e.nom_45/data_200715$e.nom_total
data_200715$g_13.1 <- data_200715$g_13/data_200715$g_val
data_200715$g_45.1 <- data_200715$g_45/data_200715$g_val
data_200715$Comparecimento1t.1 <- data_200715$Comparecimento1t/data_200715$Aptos


# Create vectors for tests

X2000_incomepercapita <- NA
X2000_doctor <- NA
X2000_idh_education <- NA
X2000_idh_income <- NA
X2000_idh_longevity <- NA
X2000_illiterate <- NA
X2000_infant <- NA
X2000_pop <- NA
X2000_poverty <- NA
p_13.1 <- NA
p_45.1 <- NA
f.nom_13.1 <- NA
f.nom_45.1 <- NA
e.nom_13.1 <- NA
e.nom_45.1 <- NA
g_13.1 <- NA
g_45.1 <- NA
Comparecimento1t.1 <- NA
age <- NA
sex <- NA
ftest.h <- NA

#Renaming variables
names(X2000_incomepercapita) <- "Income per capita"
names(X2000_doctor) <- "Doctors per thousand pop."
names(X2000_idh_education) <- "Education (IDH)"
names(X2000_idh_income) <- "Income (IDH)"
names(X2000_idh_longevity) <- "Longevity (IDH)"
names(X2000_illiterate) <- "Illiteracy rate"
names(X2000_infant) <- "Infant mortality"
names(X2000_pop)  <- "Population"
names(X2000_poverty) <- "Poverty rate"
names(p_13.1) <- "Vote for Lula in 1998"
names(p_45.1) <- "Vote for FHC in 1998"
names(f.nom_13.1) <- "Vote for PT (fed. dep.) in 1998"
names(f.nom_45.1) <- "Votes for PSDB (fed. dep.) in 1998"
names(g_13.1) <- "Vote for PT (governor) in 1998"
names(g_45.1) <- "Vote for PSDB (governor) in 1998"
names(e.nom_13.1) <- "Votes for PT (state dep.) in 1998"
names(e.nom_45.1) <- "Votes for PSBD (state dep.) in 1998"
names(Comparecimento1t.1) <- "Turnout in 1998"
names(sex) <- "Mayor's Sex"
names(age) <- "Mayor's Age"
names(ftest.h) <-  "F-test"

#Creating Standardized Variables
data_200315$X2000_incomepercapita.2 <- (data_200315$X2000_incomepercapita - mean(data_200315$X2000_incomepercapita, na.rm=T))/sd(data_200315$X2000_incomepercapita, na.rm=T)
data_200315$X2000_doctor.2 <- (data_200315$X2000_doctor- mean(data_200315$X2000_doctor, na.rm=T))/sd(data_200315$X2000_doctor, na.rm=T)
data_200315$X2000_idh_education.2 <- (data_200315$X2000_idh_education - mean(data_200315$X2000_idh_education, na.rm=T))/sd(data_200315$X2000_idh_education, na.rm=T)
data_200315$X2000_idh_income.2 <- (data_200315$X2000_idh_income - mean(data_200315$X2000_idh_income, na.rm=T))/sd(data_200315$X2000_idh_income, na.rm=T)
data_200315$X2000_idh_longevity.2 <- (data_200315$X2000_idh_longevity - mean(data_200315$X2000_idh_longevity, na.rm=T))/sd(data_200315$X2000_idh_longevity, na.rm=T)
data_200315$X2000_illiterate.2 <- (data_200315$X2000_illiterate - mean(data_200315$X2000_illiterate, na.rm=T))/sd(data_200315$X2000_illiterate, na.rm=T)
data_200315$X2000_infant.2 <- (data_200315$X2000_infant - mean(data_200315$X2000_infant, na.rm=T))/sd(data_200315$X2000_infant, na.rm=T)
data_200315$X2000_pop.2 <- (data_200315$X2000_pop - mean(data_200315$X2000_pop, na.rm=T))/sd(data_200315$X2000_pop, na.rm=T)
data_200315$X2000_poverty.2 <- (data_200315$X2000_poverty - mean(data_200315$X2000_poverty, na.rm=T))/sd(data_200315$X2000_poverty, na.rm=T)
data_200315$p_13.2 <- (data_200315$p_13.1 - mean(data_200315$p_13.1, na.rm=T))/sd(data_200315$p_13.1, na.rm=T)
data_200315$p_45.2 <- (data_200315$p_45.1 - mean(data_200315$p_45.1, na.rm=T))/sd(data_200315$p_45.1, na.rm=T)
data_200315$f.nom_13.2 <- (data_200315$f.nom_13.1 - mean(data_200315$f.nom_13.1, na.rm=T))/sd(data_200315$f.nom_13.1, na.rm=T) 
data_200315$f.nom_45.2 <- (data_200315$f.nom_45.1 - mean(data_200315$f.nom_45.1, na.rm=T))/sd(data_200315$f.nom_45.1, na.rm=T)
data_200315$g_13.2 <- (data_200315$g_13.1 - mean(data_200315$g_13.1, na.rm=T))/sd(data_200315$g_13.1, na.rm=T) 
data_200315$g_45.2 <- (data_200315$g_45.1 - mean(data_200315$g_45.1, na.rm=T))/sd(data_200315$g_45.1, na.rm=T)
data_200315$e.nom_13.2 <- (data_200315$e.nom_13.1 - mean(data_200315$e.nom_13.1, na.rm=T))/sd(data_200315$e.nom_13.1, na.rm=T)
data_200315$e.nom_45.2 <- (data_200315$e.nom_45.1 - mean(data_200315$e.nom_45.1, na.rm=T))/sd(data_200315$e.nom_45.1, na.rm=T)
data_200315$Comparecimento1t.2 <- (data_200315$Comparecimento1t.1 - mean(data_200315$Comparecimento1t.1, na.rm=T))/sd(data_200315$Comparecimento1t.1, na.rm=T)
data_200315$sex.2 <- (data_200315$sex - mean(data_200315$sex, na.rm=T))/sd(data_200315$sex, na.rm=T)
data_200315$age.2 <- (data_200315$age - mean(data_200315$age, na.rm=T))/sd(data_200315$age, na.rm=T)


#Creating Standardized Variables
data_200715$X2000_incomepercapita.2 <- (data_200715$X2000_incomepercapita - mean(data_200715$X2000_incomepercapita, na.rm=T))/sd(data_200715$X2000_incomepercapita, na.rm=T)
data_200715$X2000_doctor.2 <- (data_200715$X2000_doctor- mean(data_200715$X2000_doctor, na.rm=T))/sd(data_200715$X2000_doctor, na.rm=T)
data_200715$X2000_idh_education.2 <- (data_200715$X2000_idh_education - mean(data_200715$X2000_idh_education, na.rm=T))/sd(data_200715$X2000_idh_education, na.rm=T)
data_200715$X2000_idh_income.2 <- (data_200715$X2000_idh_income - mean(data_200715$X2000_idh_income, na.rm=T))/sd(data_200715$X2000_idh_income, na.rm=T)
data_200715$X2000_idh_longevity.2 <- (data_200715$X2000_idh_longevity - mean(data_200715$X2000_idh_longevity, na.rm=T))/sd(data_200715$X2000_idh_longevity, na.rm=T)
data_200715$X2000_illiterate.2 <- (data_200715$X2000_illiterate - mean(data_200715$X2000_illiterate, na.rm=T))/sd(data_200715$X2000_illiterate, na.rm=T)
data_200715$X2000_infant.2 <- (data_200715$X2000_infant - mean(data_200715$X2000_infant, na.rm=T))/sd(data_200715$X2000_infant, na.rm=T)
data_200715$X2000_pop.2 <- (data_200715$X2000_pop - mean(data_200715$X2000_pop, na.rm=T))/sd(data_200715$X2000_pop, na.rm=T)
data_200715$X2000_poverty.2 <- (data_200715$X2000_poverty - mean(data_200715$X2000_poverty, na.rm=T))/sd(data_200715$X2000_poverty, na.rm=T)
data_200715$p_13.2 <- (data_200715$p_13.1 - mean(data_200715$p_13.1, na.rm=T))/sd(data_200715$p_13.1, na.rm=T)
data_200715$p_45.2 <- (data_200715$p_45.1 - mean(data_200715$p_45.1, na.rm=T))/sd(data_200715$p_45.1, na.rm=T)
data_200715$f.nom_13.2 <- (data_200715$f.nom_13.1 - mean(data_200715$f.nom_13.1, na.rm=T))/sd(data_200715$f.nom_13.1, na.rm=T) 
data_200715$f.nom_45.2 <- (data_200715$f.nom_45.1 - mean(data_200715$f.nom_45.1, na.rm=T))/sd(data_200715$f.nom_45.1, na.rm=T)
data_200715$g_13.2 <- (data_200715$g_13.1 - mean(data_200715$g_13.1, na.rm=T))/sd(data_200715$g_13.1, na.rm=T) 
data_200715$g_45.2 <- (data_200715$g_45.1 - mean(data_200715$g_45.1, na.rm=T))/sd(data_200715$g_45.1, na.rm=T)
data_200715$e.nom_13.2 <- (data_200715$e.nom_13.1 - mean(data_200715$e.nom_13.1, na.rm=T))/sd(data_200715$e.nom_13.1, na.rm=T)
data_200715$e.nom_45.2 <- (data_200715$e.nom_45.1 - mean(data_200715$e.nom_45.1, na.rm=T))/sd(data_200715$e.nom_45.1, na.rm=T)
data_200715$Comparecimento1t.2 <- (data_200715$Comparecimento1t.1 - mean(data_200715$Comparecimento1t.1, na.rm=T))/sd(data_200715$Comparecimento1t.1, na.rm=T)
data_200715$sex.2 <- (data_200715$sex- mean(data_200715$sex, na.rm=T))/sd(data_200715$sex, na.rm=T)
data_200715$age.2 <- (data_200715$age - mean(data_200715$age, na.rm=T))/sd(data_200715$age, na.rm=T) 


#Pre-treatment covariates
#List of names
outcomes <-  list(
  data_200315$X2000_incomepercapita.2,
  data_200315$X2000_doctor.2,
  data_200315$X2000_idh_education.2,
  data_200315$X2000_idh_income.2,
  data_200315$X2000_idh_longevity.2,
  data_200315$X2000_illiterate.2,
  data_200315$X2000_infant.2,
  data_200315$X2000_pop.2, 
  data_200315$X2000_poverty.2,
  data_200315$p_13.2, 
  data_200315$p_45.2,
  data_200315$f.nom_13.2, 
  data_200315$f.nom_45.2,
  data_200315$g_13.2, 
  data_200315$g_45.2,
  data_200315$e.nom_13.2,
  data_200315$e.nom_45.2,
  data_200315$Comparecimento1t.2,
  data_200315$sex.2, 
  data_200315$age.2)

#Pre-treatment covariates
#List of names
outcomes.state <-  list(
  data_200715$X2000_incomepercapita.2,
  data_200715$X2000_doctor.2,
  data_200715$X2000_idh_education.2,
  data_200715$X2000_idh_income.2,
  data_200715$X2000_idh_longevity.2,
  data_200715$X2000_illiterate.2,
  data_200715$X2000_infant.2,
  data_200715$X2000_pop.2, 
  data_200715$X2000_poverty.2,
  data_200715$p_13.2, 
  data_200715$p_45.2,
  data_200715$f.nom_13.2, 
  data_200715$f.nom_45.2,
  data_200715$g_13.2, 
  data_200715$g_45.2,
  data_200715$e.nom_13.2,
  data_200715$e.nom_45.2,
  data_200715$Comparecimento1t.2,
  data_200715$sex.2, 
  data_200715$age.2)


names <- c("Income per capita", 
           "Doctors per thousand pop.",
           "Education (IDH)",
           "Income (IDH)", 
           "Longevity (IDH)", 
           "Illiteracy rate", 
           "Infant mortality", 
           "Population",
           "Poverty rate", 
           "Vote for Lula",
           "Vote for FHC", 
           "Vote for PT (fed. dep.)",
           "Votes for PSDB (fed. dep.)",
           "Vote for PT (governor)",
           "Vote for PSDB (governor)",
           "Votes for PT (state dep.)",
           "Votes for PSBD (state dep.)",
           "Turnout", 
           "Mayor's Sex",
           "Mayor's Age")

############################################### Tables Supplemental Material

#S.B.1 NSPs' Areas of Expertise (2003-2011) 

xtable(prop.table(table(ipea_orgs$IBGE_grupo))*100)
nrow(ipea_orgs)

#S.B.2 Types of Organizations (2003-2011)

xtable(prop.table(table(ipea_orgs$class_IPEA_reduzido))*100)
nrow(ipea_orgs)

#S.B.3 NSPs' Sources of Funding (2009-2015) 

table_sources <- as.matrix(apply(sources_orgs[,2:ncol(sources_orgs)], 2, mean, na.rm=T))
table_sourcesf <- matrix(NA, 5, 2)
table_sourcesf[,2] <- round(as.numeric(c(table_sources[1,1] + table_sources[2,1],
						table_sources[3,1] + table_sources[4,1],
						table_sources[5,1],
						table_sources[6,1] + table_sources[7,1],
						100 - table_sources[8,1])), 2)
table_sourcesf[,1] <- c("Own Funds (services and members' fees)", "Private donations", 
						"Public Funds", "International Funds", "Missing")
xtable(table_sourcesf)
nrow(sources_orgs)

#S.B.4 Number of organizations per municipality (2008-2015) 

esfl$ibge <- substr(esfl$IBGEcod, 1, nchar(esfl$IBGEcod)-1)
orgs.mun <- table(esfl$ibge, esfl$CD_IDENTIF_PROPONENTE)
orgs.mun <- cbind(orgs.mun, Total = rowSums(orgs.mun))
smr <- summary(orgs.mun[,ncol(orgs.mun)])
orgs <- table(orgs.mun[,ncol(orgs.mun)])
prop.table(orgs)*100
orgs[1]/sum(orgs)
xtable(as.matrix(smr))
length(unique(esfl$ibge))
length(unique(esfl$CD_IDENTIF_PROPONENTE))

#S.B.5 Share per Organization (2008-2015)

temp <- esfl %>% group_by(CD_IDENTIF_PROPONENTE, ibge) %>% summarise(soma = sum(d_empenhado)) %>% filter(soma >= 0)
tempm <- esfl %>% group_by(ibge) %>% summarise(soma_mun = sum(d_empenhado))
temp.all <- temp %>% left_join(tempm, by = "ibge")
temp.all <- temp.all %>% mutate(share = soma/soma_mun)
xtable(as.matrix(summary(temp.all$share)))
length(unique(temp.all$ibge))
length(unique(temp.all$CD_IDENTIF_PROPONENTE))

#S.B.6 Top 10 Ministries in Transfers to NSPs (2008–2015) 

spendinge <- with(esfl, aggregate(d_empenhado, by = list(NM_ORGAO_CONCEDENTE), sum, na.rm = T))
spendinge <- spendinge[order(spendinge$x, decreasing = T),]
xtable(spendinge[1:10,])

#S.B.7 Top 10 Ministries in Transfers to Mayors (2008–2015)

mayors <- sic[which(sic$TX_ESFERA_ADM_PROPONENTE == "MUNICIPAL"),]
spendingm <- with(mayors, aggregate(d_empenhado, by= list(NM_ORGAO_CONCEDENTE), sum, na.rm = T))
spendingm <- spendingm[order(spendingm$x, decreasing = T),]
xtable(spendingm[1:10,])

#S.C.1 Descriptive Statistics: Municipalities in which the aligned mayor was either a 
#winner or runner-up, 2003-2015 

d_summary <- data_200315[c("vote_margin_share2", "dtransfers.pc", 
							"dtransfers.pc.mayors",
                            "VOTO_MUN_TOTAL", "X2000_incomepercapita", 
                            "X2000_doctor", 
                            "X2000_idh_education", "X2000_idh_income",
							"X2000_idh_longevity", "X2000_illiterate", 
							"X2000_infant", "X2000_pop",
							"X2000_poverty", "p_13.1", "p_45.1",
							"f.nom_13.1", "f.nom_45.1", "g_13.1", 
							"g_45.1", "e.nom_13.1", "e.nom_45.1", 
							"Comparecimento1t.1", "sex", "age")]

xtable(summary.stats(d_summary))

#S.C.2 Balance Tests, Difference of Means

#Excluding missing data because will use randomization inference
data_200315v2 <- data_200315
data_200315v3 <- data_200315v2 %>% filter(!is.na(X2000_incomepercapita.2))
cat(nrow(data_200315v2) - nrow(data_200315v3), "cases missing: demographics") #change this back to v3
data_200315v4 <- data_200315v3 %>% filter(!is.na(age.2))
cat(nrow(data_200315v3) - nrow(data_200315v4), "cases missing: age")
rm(data_200315v2, data_200315v3)

#Getting pre-treatment outcomes
outcomes2 <-  list(
  data_200315v4$X2000_incomepercapita.2,
  data_200315v4$X2000_doctor.2,
  data_200315v4$X2000_idh_education.2,
  data_200315v4$X2000_idh_income.2,
  data_200315v4$X2000_idh_longevity.2,
  data_200315v4$X2000_illiterate.2,
  data_200315v4$X2000_infant.2,
  data_200315v4$X2000_pop.2, 
  data_200315v4$X2000_poverty.2,
  data_200315v4$p_13.2, 
  data_200315v4$p_45.2,
  data_200315v4$f.nom_13.2, 
  data_200315v4$f.nom_45.2,
  data_200315v4$g_13.2, 
  data_200315v4$g_45.2,
  data_200315v4$e.nom_13.2,
  data_200315v4$e.nom_45.2,
  data_200315v4$Comparecimento1t.2,
  data_200315v4$sex.2,
  data_200315v4$age.2)

names <- c("Income per capita", 
           "Doctors per thousand pop.",
           "Education (IDH)",
           "Income (IDH)", 
           "Longevity (IDH)", 
           "Illiteracy rate", 
           "Infant mortality", 
           "Population",
           "Poverty rate", 
           "Vote for Lula",
           "Vote for FHC", 
           "Vote for PT (fed. dep.)",
           "Votes for PSDB (fed. dep.)",
           "Vote for PT (governor)",
           "Vote for PSDB (governor)",
           "Votes for PT (state dep.)",
           "Votes for PSBD (state dep.)",
           "Turnout", 
           "Mayor's Sex",
           "Mayor's Age")

#Very slow, be aware
balance_table <- rdd.table.balance(runvar = data_200315v4$vote_margin_share, 
                                   treat = data_200315v4$treat,
                                   outcome = outcomes2, title=names) 

balance_table <- do.call(rbind, balance_table)
xtable(balance_table)

#S.C.3 Balance Tests, local polynomial regression

poly_balance_table <- rdd.table.balance.poly(runvar = data_200315$vote_margin_share, outcome = outcomes)

poly_balance_table <- do.call(rbind, poly_balance_table)

rownames(poly_balance_table) <- names
xtable(poly_balance_table[,c(1,2)])

#S.C.4 RDD Table Lagged outcomes

#Clusters
data_lagged$ibge2 <- as.numeric(data_lagged$ibge)
#Nonzero
data_lagged$nonzero <- ifelse(data_lagged$dtransfers.pc > 0, 1, 0)
data_lagged$nonzero.mayors <- ifelse(data_lagged$dtransfers.pc.mayors > 0, 1, 0)

#Optimal Bandwiths using CCT functions
nsp_cctbwl <- rdbwselect(y = data_lagged$dtransfers.pc,  x= data_lagged$vote_margin_share,
                          all = TRUE, kernel = "uniform", p=  1, q = 2)
mayors_cctbwl <- rdbwselect(y = data_lagged$dtransfers.pc.mayors, x = data_lagged$vote_margin_share, 
                            all = TRUE, kernel = "tri", p = 1, q = 2)

opt.bw.nspl <- optimal.bw(Y = data_lagged$dtransfers.pc, 
                         X = data_lagged$vote_margin_share, 
                         c = 0, reg = T)

opt.bw.mayorsl <- optimal.bw(Y = data_lagged$dtransfers.pc.mayors, 
                            X = data_lagged$vote_margin_share, 
                            c = 0, reg = T)

rdd.table(runvar=data_lagged$vote_margin_share, treat=data_lagged$treat,
          outcome.1= data_lagged$dtransfers.pc, clusterT=T, cluster.var=data_lagged$ibge2,
          non.zero= data_lagged$nonzero, p=1, q=2, kernel="uniform",
          opt.pc= opt.bw.nspl)

rdd.table(runvar=data_lagged$vote_margin_share, treat= data_lagged$treat,
          outcome.1= data_lagged$dtransfers.pc.mayors, clusterT= T, cluster.var=data_lagged$ibge2,
          non.zero= data_lagged$nonzero.mayors, p=1, q=2, kernel= "tri",
          opt.pc= opt.bw.mayorsl)
          
#S.C.5 Manipulation tetss

#McCrary
DCdensity(data_200315$vote_margin_share)

# manipulation test with jackknife standard error
margin <- data_200315$vote_margin_share
rddensity(X = margin, vce="jackknife",  p = 1, q = 2)
rddensity(X = margin, kernel = "uniform", p = 1, q = 2)
rddensity(X = margin, fitselect="restricted", vce="plugin", p = 1, q = 2)

results_manipulation <- matrix(NA, 1, 3)
results_manipulation[,1:3] <- c(0.735, 0.829, 0.6148)
colnames(results_manipulation) <- c("Robust, unrestricted, triangular",
"Robust, unrestricted, uniform", "Robust, restricted, triangular")
          
         
#S.C.6 Placebo Treatments NSP

nsp_cctbw <- rdbwselect(y = data_200315$dtransfers.pc,  x= data_200315$vote_margin_share,
                        all = TRUE, kernel = "uniform", p=  1, q = 2)

table_nsp <- matrix(NA, 4, 6)

table_nsp[1,] <- round(as.numeric(summary(with(data_200315, rdrobust(y = dtransfers.pc, x = vote_margin_share, 
				   c = 0.08555076/2, kernel = "uniform")))$coefficients[3,]), 2)
table_nsp[2,] <- round(as.numeric(summary(with(data_200315, rdrobust(y = dtransfers.pc, x = vote_margin_share, 
						  c = -0.08555076/2, kernel = "uniform")))$coefficients[3,]), 2)
table_nsp[3,] <- round(as.numeric(summary(with(data_200315, rdrobust(y = dtransfers.pc, x = vote_margin_share, 
						   c = 0.08555076*2, kernel = "uniform")))$coefficients[3,]), 2)
table_nsp[4,] <- round(as.numeric(summary(with(data_200315, rdrobust(y = dtransfers.pc, x = vote_margin_share, 
						   c = -0.08555076*2, kernel = "uniform")))$coefficients[3,]	), 2)		

colnames(table_nsp) <- c("Coeff", "Std. Err.", "z", "P>|z|", "CI Lower", "CI Upper") 
rownames(table_nsp) <- c("Pos. half bandwidth", "Neg. half bandwidth", "Pos. double bandwidth", "Neg. double bandwidth")
xtable(table_nsp)

#S.C.7 Placebo Treatments Mayors

mayors_cctbw <- rdbwselect(y = data_200315$dtransfers.pc.mayors, x = data_200315$vote_margin_share,
                           all = TRUE, kernel = "tri", p = 1, q = 2)

table_mayors <- matrix(NA, 4, 6)

table_mayors[1,] <- round(as.numeric(summary(with(data_200315, rdrobust(y = dtransfers.pc.mayors, x = vote_margin_share, 
						  c = 0.1620125/2, kernel = "tri")))$coefficients[3,]), 2)
table_mayors[2,] <- round(as.numeric(summary(with(data_200315, rdrobust(y = dtransfers.pc.mayors, x = vote_margin_share, 
						  c = -0.1620125/2, kernel = "tri")))$coefficients[3,]), 2)
table_mayors[3,] <- round(as.numeric(summary(with(data_200315, rdrobust(y = dtransfers.pc.mayors, x = vote_margin_share, 
						   c = 0.1620125*2, kernel = "tri")))$coefficients[3,]), 2)
table_mayors[4,] <- round(as.numeric(summary(with(data_200315, rdrobust(y = dtransfers.pc.mayors, x = vote_margin_share, 
						   c = -0.1620125*2, kernel = "tri")))$coefficients[3,]), 2)					   

colnames(table_mayors) <- c("Coeff", "Std. Err.", "z", "P>|z|", "CI Lower", "CI Upper") 
rownames(table_mayors) <- c("Pos. half bandwidth", "Neg. half bandwidth", "Pos. double bandwidth", "Neg. double bandwidth")

xtable(table_mayors)

#S.C.8 Sum of RDD Estimates, NSP and Mayoral Federal Transfers, 2003–2015 (reais per capita)

H1 <- c(0.005, 0.01, 0.02, 0.03, 0.04, 0.05)*100

rep <- 10000

dif_bootstrap_matrix <- matrix(NA, rep, length(H1))

set.seed(06511)
for (j in 1:length(H1)){
  temp <- data_200315[abs(data_200315$vote_margin_share2) <= H1[j],]
  
  treat <- temp[temp$treat==1,]
  control <- temp[temp$treat==0,]
  
  diff_bootstrap <- rep(NA, rep)
  for (i in 1:rep) {
    treat_bootstrap <- treat[sample(1:nrow(treat), nrow(treat), TRUE),]
    control_bootstrap <- control[sample(1:nrow(control), nrow(control), TRUE),]
    dif_nsp <- mean(treat_bootstrap$dtransfers.pc, na.rm=T) - mean(control_bootstrap$dtransfers.pc, na.rm=T)
    dif_mayor <-   mean(treat_bootstrap$dtransfers.pc.mayors, na.rm=T) - mean(control_bootstrap$dtransfers.pc.mayors, na.rm=T)
    print(i)
    diff_bootstrap[i] = dif_mayor + dif_nsp
  }
  print(j)
  dif_bootstrap_matrix[,j]  <- diff_bootstrap
  
}

sume <- apply(dif_bootstrap_matrix, 2, mean) #checking bootstrap
stde <- apply(dif_bootstrap_matrix, 2, sd) #bootstrapped standard errors
tstat <- sume/stde
pval <- 1 - pnorm(tstat)

xtable(rbind(sume, stde, tstat, pval))

#S.C.9 Political Parties in the Governing Coalition 
#no code necessary -- data provided by Joyce Luz

#S.C.10 RDD Estimates (coalition), 2003-2015 (reais per capita) 
data_200315alignv2 <- data_200315align[(!is.na(data_200315align$dtransfers.pc) | !is.na(data_200315align$dtransfers.pc.mayors)),]

opt.bw.nspc <- optimal.bw(Y = data_200315alignv2$dtransfers.pc, X = data_200315alignv2$vote_margin_share, 
                          c = 0, reg = T)

opt.bw.mayorsc <- optimal.bw(Y = data_200315alignv2$dtransfers.pc.mayors, X = data_200315alignv2$vote_margin_share, 
                             c = 0, reg = T)

data_200315alignv2$nonzero <- ifelse(data_200315alignv2$dtransfers.pc==0, 0, 1)
data_200315alignv2$nonzero.mayors <- ifelse(data_200315alignv2$dtransfers.pc.mayors==0, 0, 1)
data_200315alignv2$ibge2 <- as.numeric(data_200315alignv2$ibge)

nsp_cctbwc <- rdbwselect(y = data_200315alignv2$dtransfers.pc,  x= data_200315alignv2$vote_margin_share,
                          all = TRUE, kernel = "uniform", p=  1, q = 2)
mayors_cctbwc <- rdbwselect(y = data_200315alignv2$dtransfers.pc, x = data_200315alignv2$vote_margin_share,
                             all = TRUE, kernel = "tri", p = 1, q = 2)


rdd.table(runvar= data_200315alignv2$vote_margin_share, treat=data_200315alignv2$treat,
               outcome.1= data_200315alignv2$dtransfers.pc, clusterT=T, cluster.var=data_200315alignv2$ibge2,
               non.zero= data_200315alignv2$nonzero, p=1, q=2, kernel="uniform",
               opt.pc= opt.bw.nspc)

rdd.table(runvar= data_200315alignv2$vote_margin_share, treat=data_200315alignv2$treat,
          outcome.1= data_200315alignv2$dtransfers.pc.mayors, clusterT=T, cluster.var=data_200315alignv2$ibge2,
          non.zero= data_200315alignv2$nonzero, p=1, q=2, kernel="tri",
          opt.pc= opt.bw.nspc)

#S.C.11 RDD Estimates (ministries run by the Workers’ Party), NSP Transfers, 2003-2015 (reais per capita) 

data_200315nsp_minv2 <- data_200315nsp_min[!is.na(data_200315nsp_min$dtransfers.pc),]
data_200315nsp_minv2$nonzero <- ifelse(data_200315nsp_minv2$dtransfers.pc==0, 0, 1)
data_200315nsp_minv2$ibge2 <- as.numeric(data_200315nsp_minv2$ibge)

opt.bw.nsppt_min <- optimal.bw(Y = data_200315nsp_minv2$dtransfers.pc, X = data_200315nsp_minv2$vote_margin_share, 
                          c = 0, reg = T)

data_200915mayors_minv2 <- data_200915mayors_min[!is.na(data_200915mayors_min$dtransfers.pc),]
data_200915mayors_minv2$nonzero <- ifelse(data_200915mayors_minv2$dtransfers.pc==0, 0, 1)
data_200915mayors_minv2$ibge2 <- as.numeric(data_200915mayors_minv2$ibge)

opt.bw.mayorspt_min <- optimal.bw(Y = data_200915mayors_minv2$dtransfers.pc, X = data_200915mayors_minv2$vote_margin_share, 
                               c = 0, reg = T)

nsp_cctbwcmin <- rdbwselect(y = data_200315nsp_minv2$dtransfers.pc,  x= data_200315nsp_minv2$vote_margin_share,
                            all = TRUE, kernel = "uniform", p=  1, q = 2)
mayors_cctbwmin <- rdbwselect(y = data_200915mayors_minv2$dtransfers.pc, 
                              x = data_200915mayors_minv2$vote_margin_share,
                              all = TRUE, kernel = "tri", p = 1, q = 2)

rdd.table(runvar= data_200315nsp_minv2$vote_margin_share, treat= data_200315nsp_minv2$treat,
          outcome.1= data_200315nsp_minv2$dtransfers.pc, clusterT=T, cluster.var = data_200315nsp_minv2$ibge2,
          non.zero= data_200315nsp_minv2$nonzero, p=1, q=2, kernel="uniform",
          opt.pc= opt.bw.nsppt_min)

rdd.table(runvar= data_200915mayors_minv2$vote_margin_share, treat= data_200915mayors_minv2$treat,
          outcome.1= data_200915mayors_minv2$dtransfers.pc, clusterT=T, cluster.var = data_200915mayors_minv2$ibge2,
          non.zero= data_200915mayors_minv2$nonzero, p=1, q=2, kernel="tri",
          opt.pc= opt.bw.mayorspt_min)
          
          
#S.C.12 RDD Estimates (ministries run by coalition parties), NSP Transfers, 2009-2015 (reais per capita)

data_200315nsp_notptv2 <- data_200315nsp_notpt[!is.na(data_200315nsp_notpt$dtransfers.pc),]
data_200315nsp_notptv2$nonzero <- ifelse(data_200315nsp_notptv2$dtransfers.pc==0, 0, 1)
data_200315nsp_notptv2$ibge2 <- as.numeric(data_200315nsp_notptv2$ibge)

opt.bw.nsp_npt <- optimal.bw(Y = data_200315nsp_notptv2$dtransfers.pc, X = data_200315nsp_notptv2$vote_margin_share, 
                             c = 0, reg = T)

data_200915mayors_notptv2 <- data_200915mayors_notpt[!is.na(data_200915mayors_notpt$dtransfers.pc),]
data_200915mayors_notptv2$nonzero <- ifelse(data_200915mayors_notptv2$dtransfers.pc==0, 0, 1)
data_200915mayors_notptv2$ibge2 <- as.numeric(data_200915mayors_notptv2$ibge)

opt.bw.mayors_npt <- optimal.bw(Y = data_200915mayors_notptv2$dtransfers.pc, X = data_200915mayors_notptv2$vote_margin_share, 
                                c = 0, reg = T)

nsp_cctbwcn_min <- rdbwselect(y = data_200315nsp_notptv2$dtransfers.pc,  x= data_200315nsp_notptv2$vote_margin_share,
                              all = TRUE, kernel = "uniform", p=  1, q = 2)
mayors_cctbwn_min <- rdbwselect(y = data_200915mayors_notptv2$dtransfers.pc, 
                                x = data_200915mayors_notptv2$vote_margin_share,
                                all = TRUE, kernel = "tri", p = 1, q = 2)

rdd.table(runvar= data_200315nsp_notptv2$vote_margin_share, treat= data_200315nsp_notptv2$treat,
          outcome.1= data_200315nsp_notptv2$dtransfers.pc, clusterT=T, cluster.var = data_200315nsp_notptv2$ibge2,
          non.zero= data_200315nsp_notptv2$nonzero, p=1, q=2, kernel="uniform",
          opt.pc= opt.bw.nsp_npt)

rdd.table(runvar= data_200915mayors_notptv2$vote_margin_share, treat= data_200915mayors_notptv2$treat,
          outcome.1= data_200915mayors_notptv2$dtransfers.pc, clusterT=T, cluster.var = data_200915mayors_notptv2$ibge2,
          non.zero= data_200915mayors_notptv2$nonzero, p=1, q=2, kernel="tri",
          opt.pc= opt.bw.mayors_npt)
          

#S.C.13 RDD Estimates (below median poverty level), NSP and Mayoral Transfers, 2003-2015 (reais per capita) 

#Getting data ready
data_200315$hdi <- (data_200315$X2000_idh_education*data_200315$X2000_idh_income*data_200315$X2000_idh_longevity)^(1/3)

data_200315belowh <- data_200315[which(data_200315$hdi <= median(data_200315$hdi, na.rm=T)),]
data_200315aboveh <- data_200315[which(data_200315$hdi > median(data_200315$hdi, na.rm=T)),]
data_200315below <- data_200315[which(data_200315$X2000_poverty <= median(data_200315$X2000_poverty, na.rm=T)),]
data_200315above <- data_200315[which(data_200315$X2000_poverty > median(data_200315$X2000_poverty, na.rm=T)),]

#Getting Estimates
nsp_cctbwb <- rdbwselect(y=data_200315below$dtransfers.pc, x=data_200315below$vote_margin_share, 
                         all=TRUE, kernel="uniform", p=1, q=2)
opt.bw.nspb <- optimal.bw(Y = data_200315below$dtransfers.pc, 
                          X = data_200315below$vote_margin_share, 
                          c = 0, reg = T)
rdd.table(runvar=data_200315below$vote_margin_share, 
          treat=data_200315below$treat, clusterT = T, cluster.var = data_200315below$ibge2,
          outcome.1= data_200315below$dtransfers.pc, 
          non.zero= data_200315below$nonzero, p=1, q=2, kernel="uniform",
          opt.pc= opt.bw.nspb)
          
mayors_cctbwb <- rdbwselect(y=data_200315below$dtransfers.pc.mayors, x=data_200315below$vote_margin_share,
                            all=TRUE, kernel="tri", p=1, q=2)
opt.bw.mayorsb <- optimal.bw(Y = data_200315below$dtransfers.pc.mayors, 
                             X = data_200315below$vote_margin_share, 
                             c = 0, reg = T)

rdd.table(runvar=data_200315below$vote_margin_share, 
          treat=data_200315below$treat, clusterT = T, cluster.var = data_200315below$ibge2,
          outcome.1= data_200315below$dtransfers.pc.mayors, 
          non.zero= data_200315below$nonzero.mayors, p=1, q=2, kernel="tri",
          opt.pc= opt.bw.mayorsb)
         
          
#S.C.14 RDD Estimates (above median poverty level), NSP and Mayoral Transfers, 2003-2015 (reais per capita)  

nsp_cctbwa <- rdbwselect(y=data_200315above$dtransfers.pc, x=data_200315above$vote_margin_share, 
                         all=TRUE, kernel="uniform", p=1, q=2)
opt.bw.nspa <- optimal.bw(Y = data_200315above$dtransfers.pc, 
                          X = data_200315above$vote_margin_share, 
                          c = 0, reg = T)
rdd.table(runvar=data_200315above$vote_margin_share, 
          treat=data_200315above$treat, clusterT = T, cluster.var = data_200315above$ibge2,
          outcome.1= data_200315above$dtransfers.pc, 
          non.zero= data_200315above$nonzero, p=1, q=2, kernel="uniform",
          opt.pc= opt.bw.nspa)
          
mayors_cctbwa <- rdbwselect(y=data_200315above$dtransfers.pc.mayors, x=data_200315above$vote_margin_share,
                            all=TRUE, kernel="tri", p=1, q=2)
opt.bw.mayorsa <- optimal.bw(Y = data_200315above$dtransfers.pc.mayors, 
                             X = data_200315above$vote_margin_share, 
                             c = 0, reg = T)

rdd.table(runvar=data_200315above$vote_margin_share, 
          treat=data_200315above$treat, clusterT = T, cluster.var = data_200315above$ibge2,
          outcome.1= data_200315above$dtransfers.pc.mayors, 
          non.zero= data_200315above$nonzero.mayors, p=1, q=2, kernel="tri",
          opt.pc= opt.bw.mayorsa)          

#S.C.15 RDD Estimates (below median IDH), NSP Transfers, 2003-2015 (reais per capita) 

nsp_cctbwbh <- rdbwselect(y=data_200315belowh$dtransfers.pc, x=data_200315belowh$vote_margin_share, 
                          all=TRUE, kernel="uniform", p=1, q=2)
opt.bw.nspbh <- optimal.bw(Y = data_200315belowh$dtransfers.pc, 
                           X = data_200315belowh$vote_margin_share, 
                           c = 0, reg = T)
rdd.table(runvar=data_200315belowh$vote_margin_share, 
          treat=data_200315belowh$treat, clusterT = T, cluster.var = data_200315belowh$ibge2,
          outcome.1= data_200315belowh$dtransfers.pc, 
          non.zero= data_200315belowh$nonzero, p=1, q=2, kernel="uniform",
          opt.pc= opt.bw.nspbh)
          
mayors_cctbwbh <- rdbwselect(y=data_200315belowh$dtransfers.pc.mayors, x=data_200315belowh$vote_margin_share,
                             all=TRUE, kernel="tri", p=1, q=2)
opt.bw.mayorsbh <- optimal.bw(Y = data_200315belowh$dtransfers.pc.mayors, 
                              X = data_200315belowh$vote_margin_share, 
                              c = 0, reg = T)

rdd.table(runvar=data_200315belowh$vote_margin_share, 
          treat=data_200315belowh$treat, clusterT = T, cluster.var = data_200315above$ibge2,
          outcome.1= data_200315belowh$dtransfers.pc.mayors, 
          non.zero= data_200315belowh$nonzero.mayors, p=1, q=2, kernel="tri",
          opt.pc= opt.bw.mayorsbh)          
          

#S.C.16 RDD Estimates (above median HDI), NSP Transfers, 2003-2015 (reais per capita)

nsp_cctbwah <- rdbwselect(y=data_200315aboveh$dtransfers.pc, x=data_200315aboveh$vote_margin_share, 
                          all=TRUE, kernel="uniform", p=1, q=2)
opt.bw.nspah <- optimal.bw(Y = data_200315aboveh$dtransfers.pc, 
                           X = data_200315aboveh$vote_margin_share, 
                           c = 0, reg = T)

rdd.table(runvar=data_200315aboveh$vote_margin_share, 
          treat=data_200315aboveh$treat, clusterT = T, cluster.var = data_200315aboveh$ibge2,
          outcome.1= data_200315aboveh$dtransfers.pc, 
          non.zero= data_200315aboveh$nonzero, p=1, q=2, kernel="uniform",
          opt.pc= opt.bw.nspah)

mayors_cctbwah <- rdbwselect(y=data_200315aboveh$dtransfers.pc.mayors, x=data_200315aboveh$vote_margin_share,
                             all=TRUE, kernel="tri", p=1, q=2)
opt.bw.mayorsah <- optimal.bw(Y = data_200315aboveh$dtransfers.pc.mayors, 
                              X = data_200315aboveh$vote_margin_share, 
                              c = 0, reg = T)

rdd.table(runvar=data_200315aboveh$vote_margin_share, 
          treat=data_200315aboveh$treat, clusterT = T, cluster.var = data_200315above$ibge2,
          outcome.1= data_200315aboveh$dtransfers.pc.mayors, 
          non.zero= data_200315aboveh$nonzero, p=1, q=2, kernel="tri",
          opt.pc= opt.bw.mayorsah)

#S.C.17 RDD estimates by election cycle, NSP transfers (reais per capita) 

rdd.table.cycle(runvar=data_200315$vote_margin_share, treat=data_200315$treat,
                outcome.1= data_200315$dtransfers.pc, yeare=data_200315$ANO_ELEICAO,
                non.zero= data_200315$nonzero, p=1, q=2, kernel="uniform")

#S.C.18 RDD estimates by election cycle, mayoral transfers (reais per capita) 

rdd.table.cycle(runvar=data_200315$vote_margin_share, treat=data_200315$treat,
                outcome.1= data_200315$dtransfers.pc.mayors, yeare=data_200315$ANO_ELEICAO,
                non.zero= data_200315$nonzero.mayors, p=1, q=2, kernel="tri")


#S.C.19 RDD Estimates, NSP and Mayoral Transfers, 2003-2015 (reais per capita), without largest observation

data_200315nsp_out0 <- data_200315[-which(data_200315$dtransfers.pc==max(data_200315$dtransfers.pc, na.rm=T)),]

#Identifying municipality as:
maximum <- data_200315[which(data_200315$dtransfers.pc==max(data_200315$dtransfers.pc, na.rm =T)),]

opt.bw.nsp0 <- optimal.bw(Y = data_200315nsp_out0$dtransfers.pc, 
                         X = data_200315nsp_out0$vote_margin_share, 
                         c = 0, reg = T)

nsp_cctbw0 <- rdbwselect(y=data_200315nsp_out0$dtransfers.pc, x=data_200315nsp_out0$vote_margin_share,
                        all=TRUE, kernel="uniform", p=1, q=2)

rdd.table(runvar=data_200315nsp_out0$vote_margin_share, treat=data_200315nsp_out0$treat,
          outcome.1= data_200315nsp_out0$dtransfers.pc, clusterT = T, cluster.var = data_200315nsp_out0$ibge2,
          non.zero= data_200315nsp_out0$nonzero, p=1, q=2, kernel="uniform",
          opt.pc= opt.bw.nsp0)

data_200315mayors_out0 <- data_200315[-which(data_200315$dtransfers.pc.mayors==max(data_200315$dtransfers.pc.mayors, na.rm = T)),]

#Identifying municipality as:
maximum <- data_200315[which(data_200315$dtransfers.pc.mayors==max(data_200315$dtransfers.pc.mayors, na.rm =T)),]
#Município Santa Rosa

opt.bw.mayors0 <- optimal.bw(Y = data_200315mayors_out0$dtransfers.pc.mayors, 
                            X = data_200315mayors_out0$vote_margin_share, 
                            c = 0, reg = T)

mayors_cctbw0 <- rdbwselect(y=data_200315mayors_out0$dtransfers.pc.mayors, x=data_200315mayors_out0$vote_margin_share,
                           all=TRUE, kernel="tri", p=1, q=2)

rdd.table(runvar=data_200315mayors_out0$vote_margin_share, treat=data_200315mayors_out0$treat,
          outcome.1= data_200315mayors_out0$dtransfers.pc.mayors, clusterT = T, cluster.var = data_200315mayors_out0$ibge2,
          non.zero= data_200315mayors_out0$nonzero.mayors, p=1, q=2, kernel="triangular",
          opt.pc= opt.bw.mayors0)

#S.C.20 RDD Estimates, NSP Transfers, 2003-2015 (reais per capita), without observations larger than three standard deviations
data_200315nsp_out1 <- data_200315[which(data_200315$dtransfers.pc < 
                                           (mean(data_200315$dtransfers.pc, na.rm=T)  
                                            + sd(data_200315$dtransfers.pc, na.rm=T)*3)), ]

opt.bw.nspo1 <- optimal.bw(Y = data_200315nsp_out1$dtransfers.pc, 
                          X = data_200315nsp_out1$vote_margin_share, 
                          c = 0, reg = T)
                          
nsp_cctbw1 <- rdbwselect(y=data_200315nsp_out1$dtransfers.pc, x=data_200315nsp_out1$vote_margin_share,
                        all=TRUE, kernel="uniform", p=1, q=2)                         

rdd.table(runvar= data_200315nsp_out1$vote_margin_share, treat= data_200315nsp_out1$treat,
          outcome.1= data_200315nsp_out1$dtransfers.pc, clusterT=T, cluster.var= data_200315nsp_out1$ibge2,
          non.zero= data_200315nsp_out1$nonzero, p=1, q=2, kernel="uniform",
          opt.pc= opt.bw.nspo1)


data_200315mayors_out1 <- data_200315[which(data_200315$dtransfers.pc.mayors < 
                                              (mean(data_200315$dtransfers.pc.mayors, na.rm=T)  
                                               + sd(data_200315$dtransfers.pc.mayors, na.rm=T)*3)), ]

opt.bw.mayorso <- optimal.bw(Y = data_200315mayors_out1$dtransfers.pc.mayors, 
                             X = data_200315mayors_out1$vote_margin_share, 
                             c = 0, reg = T)
mayors_cctbw1 <- rdbwselect(y=data_200315mayors_out1$dtransfers.pc.mayors, x=data_200315mayors_out1$vote_margin_share,
                           all=TRUE, kernel="tri", p=1, q=2)                             

rdd.table(runvar= data_200315mayors_out1$vote_margin_share, treat=data_200315mayors_out1$treat,
          outcome.1= data_200315mayors_out1$dtransfers.pc.mayors, clusterT=T, cluster.var=data_200315mayors_out1$ibge2,
          non.zero= data_200315mayors_out1$nonzero.mayors, p=1, q=2, kernel="tri",
          opt.pc= opt.bw.mayorso)
          
#S.C.21 RDD Estimates, NSP Transfers, 2003-2015, logged outcomes (1+log)

#Some recoding
data_200315$ldtransfers.pc <- log(1+data_200315$dtransfers.pc)
data_200315$ldtransfers.pc.mayors <- log(1+data_200315$dtransfers.pc.mayors)

opt.bw.nspl <- optimal.bw(Y = data_200315$ldtransfers.pc, 
                          X = data_200315$vote_margin_share, 
                          c = 0, reg = T)
nsp_cctbwl <- rdbwselect(y=data_200315$ldtransfers.pc,  x=data_200315$vote_margin_share,
                             all=TRUE, kernel="uniform", p=1, q=2) 

opt.bw.mayorsl <- optimal.bw(Y = data_200315$ldtransfers.pc.mayors, 
                             X = data_200315$vote_margin_share, 
                             c = 0, reg = T)
mayors_cctbwl <- rdbwselect(y=data_200315$ldtransfers.pc.mayors,  x=data_200315$vote_margin_share,
                             all=TRUE, kernel="tri", p=1, q=2)                             

#NSP
rdd.table(runvar=data_200315$vote_margin_share, treat=data_200315$treat,
          outcome.1= data_200315$ldtransfers.pc, clusterT=T, cluster.var=data_200315$ibge2,
          non.zero= data_200315$nonzero, p=1, q=2, kernel="uniform",
          opt.pc= opt.bw.nspl)

#Mayors
rdd.table(runvar= data_200315$vote_margin_share, treat=data_200315$treat,
          outcome.1= data_200315$ldtransfers.pc.mayors, clusterT=T, cluster.var=data_200315$ibge2,
          non.zero= data_200315$nonzero.mayors, p=1, q=2, kernel="tri",
          opt.pc= opt.bw.mayorsl)


#S.C.22 RDD Estimates, NSP Transfers, 2003-2015, logged outcomes (inverse hyperbolic sine transformation)

#inverse hyperbolic sine transformation
data_200315$ldtransfers.pcinv <- inv_hyp(data_200315$dtransfers.pc)
data_200315$ldtransfers.pcinv.mayors <- inv_hyp(data_200315$dtransfers.pc.mayors)

opt.bw.nspli <- optimal.bw(Y = data_200315$ldtransfers.pcinv, 
                           X = data_200315$vote_margin_share, 
                           c = 0, reg = T)
nsp_cctbwl <- rdbwselect(y=data_200315$ldtransfers.pcinv,  x=data_200315$vote_margin_share,
                             all=TRUE, kernel="uniform", p=1, q=2) 

opt.bw.mayorsli <- optimal.bw(Y = data_200315$ldtransfers.pcinv.mayors, 
                              X = data_200315$vote_margin_share, 
                              c = 0, reg = T)
mayors_cctbwl <- rdbwselect(y=data_200315$ldtransfers.pcinv.mayors,  x=data_200315$vote_margin_share,
                             all=TRUE, kernel="tri", p=1, q=2)

#NSP
rdd.table(runvar=data_200315$vote_margin_share, treat=data_200315$treat,
          outcome.1= data_200315$ldtransfers.pcinv, clusterT=T, cluster.var=data_200315$ibge2,
          non.zero= data_200315$nonzero, p=1, q=2, kernel="uniform",
          opt.pc= opt.bw.nspli)

#Mayors
rdd.table(runvar= data_200315$vote_margin_share, treat=data_200315$treat,
          outcome.1= data_200315$ldtransfers.pcinv.mayors, clusterT=T, cluster.var=data_200315$ibge2,
          non.zero= data_200315$nonzero.mayors, p=1, q=2, kernel="tri",
          opt.pc= opt.bw.mayorsli)

#S.C.23 Descriptive Statistics: Municipalities in which the aligned mayor was either a 
#winner or runner-up, 2007-2015  (State Transfers)

d_summary <- data_200715[c("vote_margin_share2", "dtransfers.pc", 
                           "dtransfers.pc.mayors",
                           "VOTO_MUN_TOTAL", "X2000_incomepercapita", 
                           "X2000_doctor", 
                           "X2000_idh_education", "X2000_idh_income",
                           "X2000_idh_longevity", "X2000_illiterate", 
                           "X2000_infant", "X2000_pop",
                           "X2000_poverty", "p_13.1", "p_45.1",
                           "f.nom_13.1", "f.nom_45.1", "g_13.1", 
                           "g_45.1", "e.nom_13.1", "e.nom_45.1", 
                           "Comparecimento1t.1", "sex", "age")]

xtable(summary.stats(d_summary))

#S.C.24 Balance Tests, Difference of Means (State Transfers)

names <- c("Income per capita", 
           "Doctors per thousand pop.",
           "Education (IDH)",
           "Income (IDH)", 
           "Longevity (IDH)", 
           "Illiteracy rate", 
           "Infant mortality", 
           "Population",
           "Poverty rate", 
           "Vote for Lula",
           "Vote for FHC", 
           "Vote for PT (fed. dep.)",
           "Votes for PSDB (fed. dep.)",
           "Vote for PT (governor)",
           "Vote for PSDB (governor)",
           "Votes for PT (state dep.)",
           "Votes for PSBD (state dep.)",
           "Turnout", 
           "Mayor's Sex",
           "Mayor's Age")

#Be aware, it runs slow
balance_table <- rdd.table.balance(runvar = data_200715$vote_margin_share, 
                                   treat = data_200715$treat,
                                   outcome = outcomes.state, title=names) 

balance_table <- do.call(rbind, balance_table)
xtable(balance_table)

#S.C.25 Balance Tests, local polynomial regression

poly_balance_table <- rdd.table.balance.poly(runvar = data_200715$vote_margin_share, outcome = outcomes.state)

poly_balance_table <- do.call(rbind, poly_balance_table)

rownames(poly_balance_table) <- names
xtable(poly_balance_table[,c(1,2)])

#S.C.26 Tests of Manipulation of the Running Variable

DCdensity(data_200715$vote_margin_share)

margin <- data_200715$vote_margin_share
# manipulation test with jackknife standard error
rddensity(X = margin, vce="jackknife", p = 1, q = 2)
rddensity(X = margin, kernel = "uniform", p = 1, q = 2)
rddensity(X = margin, fitselect="restricted", vce="plugin", p = 1, q = 2)

results_manipulation <- matrix(NA, 1, 3)
results_manipulation[,1:3] <- c(0.5932, 0.9414, 0.0027)
colnames(results_manipulation) <- c("Robust, unrestricted, triangular",
                                    "Robust, unrestricted, uniform", "Robust, restricted, triangular")

#S.C.27 Placebo Treatments, NSP state transfers

nsp_cctbw.s <- rdbwselect(y = data_200715$dtransfers.pc,  x= data_200715$vote_margin_share,
                          all = TRUE, kernel = "uniform", p=  1, q = 2)

table_nsp.s <- matrix(NA, 4, 6)

table_nsp.s[1,] <- round(as.numeric(summary(with(data_200715, rdrobust(y = dtransfers.pc, x = vote_margin_share, 
				   c = 0.10287935/2, kernel = "uniform")))$coefficients[3,]), 2)
table_nsp.s[2,] <- round(as.numeric(summary(with(data_200715, rdrobust(y = dtransfers.pc, x = vote_margin_share, 
						  c = -0.10287935/2, kernel = "uniform")))$coefficients[3,]), 2)
table_nsp.s[3,] <- round(as.numeric(summary(with(data_200715, rdrobust(y = dtransfers.pc, x = vote_margin_share, 
						   c = 0.10287935*2, kernel = "uniform")))$coefficients[3,]), 2)
table_nsp.s[4,] <- round(as.numeric(summary(with(data_200715, rdrobust(y = dtransfers.pc, x = vote_margin_share, 
						   c = -0.10287935*2, kernel = "uniform")))$coefficients[3,]	), 2)					   

colnames(table_nsp.s) <- c("Coeff", "Std. Err.", "z", "P>|z|", "CI Lower", "CI Upper") 
rownames(table_nsp.s) <- c("Pos. half bandwidth", "Neg. half bandwidth", "Pos. double bandwidth", "Neg. double bandwidth")
xtable(table_nsp.s)

#S.C.28 Placebo Treatments, mayoral state transfers

mayors_cctbw.s <- rdbwselect(y = data_200715$dtransfers.pc.mayors, x = data_200715$vote_margin_share,
                             all = TRUE, kernel = "tri", p = 1, q = 2)

table_mayors.s <- matrix(NA, 4, 6)

table_mayors.s[1,] <- round(as.numeric(summary(with(data_200715, rdrobust(y = dtransfers.pc.mayors, x = vote_margin_share, 
						  c = 0.12660714/2, kernel = "tri")))$coefficients[3,]), 2)
table_mayors.s[2,] <- round(as.numeric(summary(with(data_200715, rdrobust(y = dtransfers.pc.mayors, x = vote_margin_share, 
						  c = -0.12660714/2, kernel = "tri")))$coefficients[3,]), 2)
table_mayors.s[3,] <- round(as.numeric(summary(with(data_200715, rdrobust(y = dtransfers.pc.mayors, x = vote_margin_share, 
						   c = 0.12660714*2, kernel = "tri")))$coefficients[3,]), 2)
table_mayors.s[4,] <- round(as.numeric(summary(with(data_200715, rdrobust(y = dtransfers.pc.mayors, x = vote_margin_share, 
						   c = -0.12660714*2, kernel = "tri")))$coefficients[3,]), 2)					   

colnames(table_mayors.s) <- c("Coeff", "Std. Err.", "z", "P>|z|", "CI Lower", "CI Upper") 
rownames(table_mayors.s) <- c("Neg. half bandwidth", "Pos. half bandwidth", "Neg. double bandwidth", "Pos. double bandwidth")
xtable(table_mayors.s)

#S.C.29 Sum of Estimates, State Transfers, 2003–2015 (reais per capita) 

H1 <- c(0.005, 0.01, 0.02, 0.03, 0.04, 0.05)*100
rep <- 10000

dif_bootstrap_matrix <- matrix(NA, rep, length(H1))

set.seed(06511)
for (j in 1:length(H1)){
  temp <- data_200715[abs(data_200715$vote_margin_share2) <= H1[j],]
  
  treat <- temp[temp$treat==1,]
  control <- temp[temp$treat==0,]
  
  diff_bootstrap <- rep(NA, rep)
  for (i in 1:rep) {
    treat_bootstrap <- treat[sample(1:nrow(treat), nrow(treat), TRUE),]
    control_bootstrap <- control[sample(1:nrow(control), nrow(control), TRUE),]
    dif_nsp <- mean(treat_bootstrap$dtransfers.pc, na.rm=T) - mean(control_bootstrap$dtransfers.pc, na.rm=T)
    dif_mayor <-   mean(treat_bootstrap$dtransfers.pc.mayors, na.rm=T) - mean(control_bootstrap$dtransfers.pc.mayors, na.rm=T)
    print(i)
    diff_bootstrap[i] = dif_mayor + dif_nsp
  }
  print(j)
  dif_bootstrap_matrix[,j]  <- diff_bootstrap
  
}

sume <- apply(dif_bootstrap_matrix, 2, mean) #checking bootstrap
stde <- apply(dif_bootstrap_matrix, 2, sd) #bootstrapped standard errors 
tstat <- sume/stde
pval <- 1 - pnorm(tstat)

xtable(rbind(sume, stde, tstat, pval))

#S.C.30 RDD Estimates, State NSP and Mayoral Transfers, 2007-2015 (reais per capita), without largest observation

#Removing extremes
data_200715nsp_out0 <- data_200715[-which(data_200715$dtransfers.pc==max(data_200715$dtransfers.pc, na.rm=T)),]

#Identifying municipality as:
maximum <- data_200715[which(data_200715$dtransfers.pc==max(data_200715$dtransfers.pc, na.rm =T)),]

opt.bw.nsp0 <- optimal.bw(Y = data_200715nsp_out0$dtransfers.pc, 
                          X = data_200715nsp_out0$vote_margin_share, 
                          c = 0, reg = T)

nsp_cctbw0 <- rdbwselect(y=data_200715nsp_out0$dtransfers.pc, x=data_200715nsp_out0$vote_margin_share,
                         all=TRUE, kernel="uniform", p=1, q=2)

rdd.table(runvar=data_200715nsp_out0$vote_margin_share, treat=data_200715nsp_out0$treat,
          outcome.1= data_200715nsp_out0$dtransfers.pc, clusterT = T, cluster.var = data_200715nsp_out0$ibge2,
          non.zero= data_200715nsp_out0$nonzero, p=1, q=2, kernel="uniform",
          opt.pc= opt.bw.nsp0)

data_200715mayors_out0 <- data_200715[-which(data_200715$dtransfers.pc.mayors==max(data_200715$dtransfers.pc.mayors, na.rm = T)),]

#Identifying municipality as:
maximum <- data_200715[which(data_200715$dtransfers.pc.mayors==max(data_200715$dtransfers.pc.mayors, na.rm =T)),]

opt.bw.mayors0 <- optimal.bw(Y = data_200715mayors_out0$dtransfers.pc.mayors, 
                             X = data_200715mayors_out0$vote_margin_share, 
                             c = 0, reg = T)

mayors_cctbw0 <- rdbwselect(y=data_200715mayors_out0$dtransfers.pc.mayors, 
                            x=data_200715mayors_out0$vote_margin_share,
                            all=TRUE, kernel="tri", p=1, q=2)

rdd.table(runvar=data_200715mayors_out0$vote_margin_share, treat=data_200715mayors_out0$treat,
          outcome.1= data_200715mayors_out0$dtransfers.pc.mayors, clusterT = T, cluster.var = data_200715mayors_out0$ibge2,
          non.zero= data_200715mayors_out0$nonzero.mayors, p=1, q=2, kernel="triangular",
          opt.pc= opt.bw.mayors0)

#S.C.31 RDD Estimates, State NSP and Mayoral Transfers, 2007-2015 (reais per capita), 
# removing observations larger than three standard deviations above the mean

data_200715nsp_out1 <- data_200715[which(data_200715$dtransfers.pc < 
                                           (mean(data_200715$dtransfers.pc, na.rm=T)  
                                            + sd(data_200715$dtransfers.pc, na.rm=T)*3)), ]

opt.bw.nspo1 <- optimal.bw(Y = data_200715nsp_out1$dtransfers.pc, 
                           X = data_200715nsp_out1$vote_margin_share, 
                           c = 0, reg = T)

nsp_cctbw1 <- rdbwselect(y=data_200715nsp_out1$dtransfers.pc, x=data_200715nsp_out1$vote_margin_share,
                         all=TRUE, kernel="uniform", p=1, q=2)                         

rdd.table(runvar= data_200715nsp_out1$vote_margin_share, treat= data_200715nsp_out1$treat,
          outcome.1= data_200715nsp_out1$dtransfers.pc, clusterT=T, cluster.var= data_200715nsp_out1$ibge2,
          non.zero= data_200715nsp_out1$nonzero, p=1, q=2, kernel="uniform",
          opt.pc= opt.bw.nspo1)


data_200715mayors_out1 <- data_200715[which(data_200715$dtransfers.pc.mayors < 
                                              (mean(data_200715$dtransfers.pc.mayors, na.rm=T)  
                                               + sd(data_200715$dtransfers.pc.mayors, na.rm=T)*3)), ]

opt.bw.mayorso <- optimal.bw(Y = data_200715mayors_out1$dtransfers.pc.mayors, 
                             X = data_200715mayors_out1$vote_margin_share, 
                             c = 0, reg = T)

mayors_cctbw1 <- rdbwselect(y=data_200715mayors_out1$dtransfers.pc.mayors, x=data_200715mayors_out1$vote_margin_share,
                            all=TRUE, kernel="tri", p=1, q=2)                             

rdd.table(runvar= data_200715mayors_out1$vote_margin_share, treat=data_200715mayors_out1$treat,
          outcome.1= data_200715mayors_out1$dtransfers.pc.mayors, clusterT=T, cluster.var=data_200715mayors_out1$ibge2,
          non.zero= data_200715mayors_out1$nonzero.mayors, p=1, q=2, kernel="tri",
          opt.pc= opt.bw.mayorso)

#S.C.32 RDD Estimates, NSP Transfers, 2007-2015, logged outcomes (1+transfers)

#some recoding
data_200715$ldtransfers.pc <- log(1+data_200715$dtransfers.pc)
data_200715$ldtransfers.pc.mayors <- log(1+data_200715$dtransfers.pc.mayors)

opt.bw.nspl <- optimal.bw(Y = data_200715$ldtransfers.pc, 
                          X = data_200715$vote_margin_share, 
                          c = 0, reg = T)
nsp_cctbwl <- rdbwselect(y=data_200715$ldtransfers.pc,  
						x=data_200715$vote_margin_share,
                             all=TRUE, kernel="uniform", p=1, q=2) 

opt.bw.mayorsl <- optimal.bw(Y = data_200715$ldtransfers.pc.mayors, 
                             X = data_200715$vote_margin_share, 
                             c = 0, reg = T)
mayors_cctbwl <- rdbwselect(y=data_200715$ldtransfers.pc.mayors,  	 x=data_200715$vote_margin_share,
                             all=TRUE, kernel="tri", p=1, q=2)                             

#NSP
rdd.table(runvar=data_200715$vote_margin_share, treat=data_200715$treat,
          outcome.1= data_200715$ldtransfers.pc, clusterT=T, cluster.var=data_200715$ibge2,
          non.zero= data_200715$nonzero, p=1, q=2, kernel="uniform",
          opt.pc= opt.bw.nspl)

#Mayors
rdd.table(runvar= data_200715$vote_margin_share, treat=data_200715$treat,
          outcome.1= data_200715$ldtransfers.pc.mayors, clusterT=T, cluster.var=data_200715$ibge2,
          non.zero= data_200715$nonzero.mayors, p=1, q=2, kernel="tri",
          opt.pc= opt.bw.mayorsl)

#S.C.33 RDD Estimates, NSP Transfers, 2007-2015, logged outcomes (inverse hyperbolic sine transformation)

#some recoding, inverse hyperbolic sine transformation
data_200715$ldtransfers.pcinv <- inv_hyp(data_200715$dtransfers.pc)
data_200715$ldtransfers.pcinv.mayors <- inv_hyp(data_200715$dtransfers.pc.mayors)

opt.bw.nspli <- optimal.bw(Y = data_200715$ldtransfers.pcinv, 
                           X = data_200715$vote_margin_share, 
                           c = 0, reg = T)
nsp_cctbwl <- rdbwselect(y=data_200715$ldtransfers.pcinv,  x=data_200715$vote_margin_share,
                             all=TRUE, kernel="uniform", p=1, q=2) 

opt.bw.mayorsli <- optimal.bw(Y = data_200715$ldtransfers.pcinv.mayors, 
                              X = data_200715$vote_margin_share, 
                              c = 0, reg = T)
mayors_cctbwl <- rdbwselect(y=data_200715$ldtransfers.pcinv.mayors,  x=data_200715$vote_margin_share,
                             all=TRUE, kernel="tri", p=1, q=2)

#NSP
rdd.table(runvar=data_200715$vote_margin_share, treat=data_200715$treat,
          outcome.1= data_200715$ldtransfers.pcinv, clusterT=T, cluster.var=data_200715$ibge2,
          non.zero= data_200715$nonzero, p=1, q=2, kernel="uniform",
          opt.pc= opt.bw.nspli)

#Mayors
rdd.table(runvar= data_200715$vote_margin_share, treat=data_200715$treat,
          outcome.1= data_200715$ldtransfers.pcinv.mayors, clusterT=T, cluster.var=data_200715$ibge2,
          non.zero= data_200715$nonzero.mayors, p=1, q=2, kernel="tri",
          opt.pc= opt.bw.mayorsli)

#S.D.1 Survey Sample Descriptive Information – Sex 

#excluding subjects less than 18 years old
pnad1 <- subset(pnad, V8005 >=18)

pnad.gender <- prop.table(table(pnad1$V0302))*100
prop.gender <- prop.table(table(obsff$gender))*100
table.gender <- cbind(prop.gender, c(pnad.gender[2], pnad.gender[1]))
table.gender <- rbind(table.gender,  c(length(!is.na(obsff$gender)), length(!is.na(pnad1$V0302))))
xtable(table.gender)

#S.D.2 Survey Sample Descriptive Information – Schooling 

pnad1$educ[pnad1$V4745==1 | pnad1$V4745==2] <- "Incomplete High School Education"
pnad1$educ[pnad1$V4745==3 | pnad1$V4745==4] <- "Incomplete High School Education"
pnad1$educ[pnad1$V4745==5 | pnad1$V4745==6] <- "Incomplete Higher Education"
pnad1$educ[pnad1$V4745==7 | pnad1$V4745==8] <- "Complete Higher Education"

obsff$educationf[obsff$education=="Analfabeto ou até a terceira série do ensino fundamental (primário incompleto)"] <- "Incomplete High School"
obsff$educationf[obsff$education=="Ensino fundamental completo e médio incompleto (ginasial completo e colegial incompleto)"] <- "Incomplete High School"
obsff$educationf[obsff$education=="Ensino médio completo e superior incompleto (colegial completo e superior incompleto)"] <- "Incomplete College"
obsff$educationf[obsff$education=="Ensino superior completo"] <- "College Degree"
obsff$educationf[obsff$education=="Pós-graduação, mestrado ou doutorado"] <- "College Degree"
obsff$educationf[obsff$education=="Quarta a sétima série do ensino fundamental (primário completo e ginásio incompleto)"] <- "Incomplete High School"

pnad.educ <- prop.table(table(pnad1$educ))*100
prop.educ <- prop.table(table(obsff$educationf))*100
table.educ <- rbind(prop.educ, pnad.educ)
table.educ <-  cbind(table.educ,  c(length(!is.na(obsff$educationf)), length(!is.na(pnad1$educ))))
xtable(t(table.educ))

#S.D.3 Survey Sample Descriptive Information – Age 

pnad1$ageg[pnad1$V8005 >= 18 & pnad1$V8005 <= 24] <- "18 to 24 years"
pnad1$ageg[pnad1$V8005 >= 25 & pnad1$V8005 <= 34]  <- "25 to 34 years"
pnad1$ageg[pnad1$V8005 >= 35 & pnad1$V8005 <= 44]  <- "35 to 44 years"
pnad1$ageg[pnad1$V8005 >= 45 & pnad1$V8005 <= 54]  <- "45 to 54 years"
pnad1$ageg[pnad1$V8005 >= 55 & pnad1$V8005 <= 64]  <- "55 to 64 years"
pnad1$ageg[pnad1$V8005 >= 65]  <- "65 years or older"

obsff$agef <- 2015 - as.numeric(obsff$age)
obsff$ageg[obsff$agef >= 18 & obsff$agef <= 24] <- "18 to 24 years"
obsff$ageg[obsff$agef >= 25 & obsff$agef <= 34]  <- "25 to 34 years"
obsff$ageg[obsff$agef >= 35 & obsff$agef <= 44]  <- "35 to 44 years"
obsff$ageg[obsff$agef >= 45 & obsff$agef <= 54]  <- "45 to 54 years"
obsff$ageg[obsff$agef >= 55 & obsff$agef <= 64]  <- "55 to 64 years"
obsff$ageg[obsff$agef >= 65]  <- "65 years or older"

pnad.age <- prop.table(table(pnad1$ageg))*100
prop.age <- prop.table(table(obsff$ageg))*100
age.table <- cbind(prop.age, pnad.age)
age.table <- rbind(age.table,  c(length(!is.na(obsff$ageg)), length(!is.na(pnad1$ageg))))
xtable(age.table)

#S.D.4 Balance Tests 

conjointff$partyf[conjointff$party=="PSDB (Partido da Social Demoracia Brasileira)"] <- 1
conjointff$partyf[conjointff$party=="PT (Partido dos Trabalhadores)"] <- 0

conjointff <- conjointff  %>% mutate(nsp_cityff = as.factor(nsp_cityf),
							  federalff =  as.factor(federalf), partisanshipff = as.factor(partisanshipf))

b.gender <- amce(genderf ~ nsp_cityff + federalff + partisanshipff,
                 data=conjointff,
                 cluster=TRUE, respondent.id="respID")

b.MW <- amce(MWf ~ nsp_cityff + federalff + partisanshipff,
             data=conjointff,
             cluster=TRUE, respondent.id="respID")

b.party <- amce(partyf ~ nsp_cityff + federalff + partisanshipff,
             data=conjointff,
             cluster=TRUE, respondent.id="respID")

 
table <- matrix(NA, 16, 4)
  
table[,1] <- c("City Govt.", "(baseline category)", "City Govt.", "Fed. Funds",
				"NSP", "", "NSP", "Fed Funds", "No Joint Provision", "(baseline category)",
				"Yes Joint Provision", "", "Partisanship (PT)", "(baseline category)", 
				"Partisanship (PSDB)", "")
				
table[,2] <- round(c(0, 0, summary(b.gender)$amce[2, 3], as.numeric(summary(b.gender)$amce[2,4]),
			   summary(b.gender)$amce[3, 3], as.numeric(summary(b.gender)$amce[3,4]), 
  			   summary(b.gender)$amce[4, 3], as.numeric(summary(b.gender)$amce[4,4]),
  			   0, 0, summary(b.gender)$amce[1, 3], as.numeric(summary(b.gender)$amce[1,4]),
  			   0, 0, summary(b.gender)$amce[5, 3], as.numeric(summary(b.gender)$amce[5,4])), 3)
  			   
table[,3] <- round(c(0, 0, summary(b.MW)$amce[2, 3], as.numeric(summary(b.MW)$amce[2,4]),
			   summary(b.MW)$amce[3, 3], as.numeric(summary(b.MW)$amce[3,4]), 
  			   summary(b.MW)$amce[4, 3], as.numeric(summary(b.MW)$amce[4,4]),
  			   0, 0, summary(b.MW)$amce[1, 3], as.numeric(summary(b.MW)$amce[1,4]),
  			   0, 0, summary(b.MW)$amce[5, 3], as.numeric(summary(b.MW)$amce[5,4])), 3) 			   

table[,4] <- round(c(0, 0,  summary(b.party)$amce[2, 3], as.numeric(summary(b.party)$amce[2,4]),
			   summary(b.party)$amce[3, 3], as.numeric(summary(b.party)$amce[3,4]), 
  			   summary(b.party)$amce[4, 3], as.numeric(summary(b.party)$amce[4,4]),
  			   0, 0, summary(b.party)$amce[1, 3], as.numeric(summary(b.party)$amce[1,4]),
  			   0, 0, summary(b.party)$amce[5, 3], as.numeric(summary(b.party)$amce[5,4])), 3)  	
  			   
xtable(table, digits=2)

#S.D.5 Effects of Welfare and Joint Provision on Probability of Mayor being Preferred
#and Change of Mayoral Ratings

conjointff <- conjointff  %>% mutate(nsp_cityff = as.factor(nsp_cityf),
							  federalff =  as.factor(federalf), partisanshipff = as.factor(partisanshipf))
             
#Analysis by aligned and unaligned
conjointff.psdb <- conjointff %>% filter(partisanship == "O prefeito pertence ao PSDB (número 45).")
conjointff.pt <- conjointff %>% filter(partisanship == "O prefeito pertence ao  PT (número 13).")

#Choice
results.choice.psdb <- amce(d_choice ~ nsp_cityff + federalff,
                            data=conjointff.psdb,
                            cluster=TRUE, respondent.id="respID")
results.choice.pt <- amce(d_choice ~ nsp_cityff + federalff,
                          data=conjointff.pt,
                          cluster=TRUE, respondent.id="respID")

#Rate
results.rate.psdb <- amce(ratef ~ nsp_cityff + federalff,
                          data=conjointff.psdb,
                          cluster=TRUE, respondent.id="respID")
results.rate.pt <- amce(ratef ~ nsp_cityff + federalff,
                        data=conjointff.pt,
                        cluster=TRUE, respondent.id="respID")
 
table <- matrix(NA, 12, 5)
  
table[,1] <- c("City Govt.", "(baseline category)", "City Govt.", "Fed. Funds",
				"NSP", "", "NSP", "Fed Funds", "No Joint Provision", "(baseline category)",
				"Yes Joint Provision", "")
				
table[,2] <- round(c(0, 0, summary(results.choice.psdb )$amce[2, 3], as.numeric(summary(results.choice.psdb)$amce[2,4]),
			   summary(results.choice.psdb )$amce[3, 3], as.numeric(summary(results.choice.psdb)$amce[3,4]), 
  			   summary(results.choice.psdb )$amce[4, 3], as.numeric(summary(results.choice.psdb)$amce[4,4]),
  			   0, 0, summary(results.choice.psdb )$amce[1, 3], as.numeric(summary(results.choice.psdb)$amce[1,4])), 3)
  			   
table[,3] <- round(c(0, 0, summary(results.choice.pt)$amce[2, 3], as.numeric(summary(results.choice.pt)$amce[2,4]),
			   summary(results.choice.pt)$amce[3, 3], as.numeric(summary(results.choice.pt)$amce[3,4]), 
  			   summary(results.choice.pt)$amce[4, 3], as.numeric(summary(results.choice.pt)$amce[4,4]),
  			   0, 0, summary(results.choice.pt)$amce[1, 3], as.numeric(summary(results.choice.pt)$amce[1,4])), 3) 			   

table[,4] <- round(c(0, 0,  summary(results.rate.psdb)$amce[2, 3], as.numeric(summary(results.rate.psdb)$amce[2,4]),
			   summary(results.rate.psdb)$amce[3, 3], as.numeric(summary(results.rate.psdb)$amce[3,4]), 
  			   summary(results.rate.psdb)$amce[4, 3], as.numeric(summary(results.rate.psdb)$amce[4,4]),
  			   0, 0, summary(results.rate.psdb)$amce[1, 3], as.numeric(summary(results.rate.psdb)$amce[1,4])), 3)  
  			   
table[,5] <- round(c(0, 0,  summary(results.rate.pt)$amce[2, 3], as.numeric(summary(results.rate.pt)$amce[2,4]),
			   summary(results.rate.pt)$amce[3, 3], as.numeric(summary(results.rate.pt)$amce[3,4]), 
  			   summary(results.rate.pt)$amce[4, 3], as.numeric(summary(results.rate.pt)$amce[4,4]),
  			   0, 0, summary(results.rate.pt)$amce[1, 3], as.numeric(summary(results.rate.pt)$amce[1,4])), 3)  			   
  			   	
  			   
xtable(table, digits=2)
length(!is.na(conjointff.psdb$d_choice))
length(!is.na(conjointff.psdb$ratef))
length(!is.na(conjointff.pt$d_choice))
length(!is.na(conjointff.pt$ratef))
length(unique(conjointff.pt$respID))
length(unique(conjointff.psdb$respID))


#S.D.6 Responsible for Santa Casa

xtable(prop.table(table(obsff$check_sc))*100)
length(!is.na(obsff$check_sc))

#S.D.7 Responsible for APAE 

xtable(prop.table(table(obsff$check_APAE))*100)
length(!is.na(obsff$check_APAE))

#S.D.8 Responsible for Municipal Hospital

xtable(prop.table(table(obsff$check_pref))*100)
length(!is.na(obsff$check_pref))

#S.D.9 Funds for Religious Hospitals

xtable(prop.table(table(obsff$Fund_hosp))*100)
length(!is.na(obsff$Fund_hosp))

#S.D.10 Funds for Municipal Hospital 

xtable(prop.table(table(obsff$Fund_local))*100)
length(!is.na(obsff$Fund_local))

#S.D.11 Funds for Welfare Associations 

xtable(prop.table(table(obsff$Fund_APAE))*100)
length(!is.na(obsff$Fund_APAE))

#S.D.12 Vignette Design
#No code necessary (see text)

#S.D.13 Balance Test – Education
vignetteff$educationf[vignetteff$education=="Analfabeto ou até a terceira série do ensino fundamental (primário incompleto)"] <- "Incomplete Primary"
vignetteff$educationf[vignetteff$education=="Ensino fundamental completo e médio incompleto (ginasial completo e colegial incompleto)"] <- "Incomplete High School"
vignetteff$educationf[vignetteff$education=="Ensino médio completo e superior incompleto (colegial completo e superior incompleto)"] <- "Incomplete College"
vignetteff$educationf[vignetteff$education=="Ensino superior completo"] <- "College Degree"
vignetteff$educationf[vignetteff$education=="Pós-graduação, mestrado ou doutorado"] <- "Graduate Degree"
vignetteff$educationf[vignetteff$education=="Quarta a sétima série do ensino fundamental (primário completo e ginásio incompleto)"] <- "Incomplete High School"

table.educ <- table(vignetteff$treat_c, vignetteff$educationf)
chisq.test(table.educ, correct = TRUE)
xtable(prop.table(table.educ, 2)*100)
table(vignetteff$educationf)

#S.D.14 Balance Test – Vote in 2014 Elections

table.vote <- table(vignetteff$treat_c, vignetteff$vote2014)
chisq.test(table.vote, correct = TRUE)
xtable(prop.table(table.vote, 2)*100)
table(vignetteff$vote2014)

#S.D.15 Balance Test – Sex

table.gender <- table(vignetteff$treat_c, vignetteff$gender)
chisq.test(table.gender, correct = TRUE)
xtable(prop.table(table.gender, 2)*100)
table(vignetteff$gender)

#S.D.16 Effects of Type of Provider and Source of Funds on Attribution of Responsibility

vignetteff$treat_prov_fed[vignetteff$treat_b== "treatb_nsp" | vignetteff$treat_b== "treatb_pref" |
                          vignetteff$treat_c== "treatc_nsp" | vignetteff$treat_c== "treatc_pref"] <- 1
vignetteff$treat_prov_fed[vignetteff$treat_b== "treatb_nsp2" | vignetteff$treat_b== "treatb_pref2" |
                          vignetteff$treat_c== "treatc_nsp2" | vignetteff$treat_c== "treatc_pref2" ] <- 0

#Credit 
#Recode outcome_credit
vignetteff$outcome_credit[vignetteff$outcome_credit== "O diretor da Associação Social Santa Maria" | 
                          vignetteff$outcome_credit== "O diretor da casa da Associação Social Santa Maria Clara" |
                          vignetteff$outcome_credit== "O diretor da casa de cuidado ao idoso"] <- "Director"


prop.table(table(vignetteff$treat_c, vignetteff$outcome_credit), 2)*100
tab.credit <- table(vignetteff$treat_c, vignetteff$outcome_credit)
chisq.test(tab.credit, correct = TRUE)
xtable(prop.table(tab.credit, 1)*100)
colSums(tab.credit)
  
#S.D.17 Effects of Type of Provider and Source of Funds on Mayoral Evaluation 

table <- matrix(NA, 3, 3)

#First cell
mean_ate1 <- round(mean(vignetteff[vignetteff$treat_c=="treatc_pref",]$outcome_creditrf),2)
se_ate1 <- round(std.error(vignetteff[vignetteff$treat_c=="treatc_pref",]$outcome_creditrf),2)
table[1, 1] <- c(paste(mean_ate1, se_ate1, sep=","))

#Second cell
mean_ate2 <- round(mean(vignetteff[vignetteff$treat_c=="treatc_pref2",]$outcome_creditrf),2)
se_ate2 <- round(std.error(vignetteff[vignetteff$treat_c=="treatc_pref2",]$outcome_creditrf),2)
table[1, 2] <- c(paste(mean_ate2, se_ate2, sep=","))

#Third cell
mean_ate3 <- round(mean(vignetteff[vignetteff$treat_c=="treatc_nsp",]$outcome_creditrf),2)
se_ate3 <- round(std.error(vignetteff[vignetteff$treat_c=="treatc_nsp",]$outcome_creditrf),2)
table[2, 1] <- c(paste(mean_ate3, se_ate3, sep=","))

#Fourth cell
mean_ate4 <- round(mean(vignetteff[vignetteff$treat_c=="treatc_nsp2",]$outcome_creditrf),2)
se_ate4 <- round(std.error(vignetteff[vignetteff$treat_c=="treatc_nsp2",]$outcome_creditrf),2)
table[2, 2] <- c(paste(mean_ate4, se_ate4, sep=","))

#Bottom cells (differences across type of provision)

#first cell 
mean_ate <- mean_ate1 - mean_ate3

#Standard error of difference of means
var_t <- var(vignetteff[vignetteff$treat_c=="treatc_nsp",]$outcome_creditrf)/nrow(vignetteff[vignetteff$treat_c=="treatc_nsp",])
var_c <- var(vignetteff[vignetteff$treat_c=="treatc_pref",]$outcome_creditrf)/nrow(vignetteff[vignetteff$treat_c=="treatc_pref",])
  
se_diff <- sqrt(var_t + var_c)

table[3, 1] <- c(paste(round(mean_ate, 2), round(se_diff, 2), sep=","))

#second cell
mean_ate <- mean_ate2 - mean_ate4

#Standard error of difference of means
var_t <- var(vignetteff[vignetteff$treat_c=="treatc_nsp2",]$outcome_creditrf)/nrow(vignetteff[vignetteff$treat_c=="treatc_nsp2",])
var_c <- var(vignetteff[vignetteff$treat_c=="treatc_pref2",]$outcome_creditrf)/nrow(vignetteff[vignetteff$treat_c=="treatc_pref2",])

se_diff <- sqrt(var_t + var_c)

table[3, 2] <- c(paste(round(mean_ate, 2), round(se_diff, 2), sep=","))

#Right-side cells (differences across class groups)

#first cell 
mean_ate <- mean_ate1 - mean_ate2
  
#Standard error difference of means
var_t <- var(vignetteff[vignetteff$treat_c=="treatc_nsp",]$outcome_creditrf)/nrow(vignetteff[vignetteff$treat_c=="treatc_nsp",])
var_c <- var(vignetteff[vignetteff$treat_c=="treatc_nsp2",]$outcome_creditrf)/nrow(vignetteff[vignetteff$treat_c=="treatc_nsp2",])

se_diff <- sqrt(var_t + var_c)

table[1, 3] <- c(paste(round(mean_ate, 2), round(se_diff, 2), sep=","))

#second cell
mean_ate <- mean_ate3 - mean_ate4
  
#Standard-error of difference of means
var_t <- var(vignetteff[vignetteff$treat_c=="treatc_pref",]$outcome_creditrf)/nrow(vignetteff[vignetteff$treat_c=="treatc_pref",])
var_c <- var(vignetteff[vignetteff$treat_c=="treatc_pref2",]$outcome_creditrf)/nrow(vignetteff[vignetteff$treat_c=="treatc_pref2",])

se_diff <- sqrt(var_t + var_c)
 
table[2, 3] <- c(paste(round(mean_ate, 2), round(se_diff, 2), sep=","))

xtable(table)

############################################### Figures Supplemental Material

#S.B.1 Growth of NSP transfers compared to overall budget growth, 
#2002–2010 in millions of reais (July 2012 values) 

pdf(paste0(dir, "paper_drafts/appendix/budgets_growth.pdf"), height=10, width=10)
par(mfrow=c(3,1), oma=c(0,2,0,0))
years <- c(2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010)
plot(growth[,6]$budget_fed_all, xlab="Years", ylab="", pch=20, lty=1, axes=F,
     ylim=c(0, 5), xlim = c(0, 9), col="red",
     cex.lab=1.5, type = "o") 
axis(1, at=seq(1,9, by=1), label=years)
axis(2, at=c(0, 2.5, 5), label=c("0", "2.5", "4.5"), las=2, cex.axis=1.5)
lines(growth[,7]$budget_fed_nsp, pch=0, lty=3, col="blue", cex.axis=1.5, type = "o")
legend("topleft", c("NSP transfers","Federal budget"), cex=1.5, 
       pch=c(16, 0), lty=c(1,3), col=c("red", "blue"), bty="n", y.intersp=0.9)

years <- c(2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010)
plot(growth[,4]$budget_state_all, xlab="Years", ylab="",
     pch=20, lty=1, axes=F, ylim=c(0, 5), xlim = c(0, 9), 
     col="red", cex.lab=1.5, type = "o")
axis(1, at=seq(1,9, by=1), label=years, cex.axis=1.5)
axis(2, at=c(0, 2.5, 4.5), label=c("0", "2.5", "4.5"), las=2, cex.axis=1.5)
lines(growth[,5]$budget_state_nsp, pch=0, lty=3, col="blue", type = "o")
legend("topleft", c("NSP transfers","State budgets"), cex=1.5, 
       pch=c(20, 0), lty=c(1,3), col=c("red", "blue"), bty="n", y.intersp=0.9)

years <- c(2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010)
plot(growth[,2]$budget_mun_all, xlab="Years", ylab="", xlim = c(0, 9),
      pch=16, axes=F, ylim=c(0, 4.5), col="red", cex.lab=1.5, 
     type = "o")
axis(1, at=seq(1,9, by=1), label=years, cex.axis=1.5)
axis(2, at=c(0, 2.5, 4.5), label=c("0", "2.5", "4.5"), las=2, cex.axis=1.5)
lines(growth[,3]$budget_mun_nsp, pch=0, lty=3, col="blue", type = "o")
legend("topleft", c("NSP transfers","Municipal budgets"), cex=1.5, 
       pch=c(16, 0), lty=c(1,3), col=c("red", "blue"), bty="n", y.intersp=0.9)
mtext("Growth rate (reference: 2002)", 2, outer=T, cex=1.2)
dev.off()

#S.C.1 P-value plot, difference of means 

X2000_incomepercapita <- NA
X2000_doctor <- NA
X2000_idh_education <- NA
X2000_idh_income <- NA
X2000_idh_longevity <- NA
X2000_illiterate <- NA
X2000_infant <- NA
X2000_pop <- NA
X2000_poverty <- NA
p_13.1 <- NA
p_45.1 <- NA
f.nom_13.1 <- NA
f.nom_45.1 <- NA
e.nom_13.1 <- NA
e.nom_45.1 <- NA
g_13.1 <- NA
g_45.1 <- NA
Comparecimento1t.1 <- NA
age <- NA
sex <- NA

#Figure 1. Regular pvalues
for (i in 1:length(H)){
  ## loop that first subsets our dataset 
  
  sel <-  abs(data_200315$vote_margin_share) <= H[i]
  data <- data_200315[sel,]  
  
  #getting the p-values
  X2000_incomepercapita[i] <- with(data, t.test.p(X2000_incomepercapita.2, treat))
  X2000_doctor[i] <- with(data, t.test.p(X2000_doctor.2, treat))
  X2000_idh_education[i] <- with(data, t.test.p(X2000_idh_education.2, treat))
  X2000_idh_income[i] <- with(data, t.test.p(X2000_idh_income.2, treat))
  X2000_idh_longevity[i] <- with(data, t.test.p(X2000_idh_longevity.2, treat))
  X2000_illiterate[i] <- with(data, t.test.p(X2000_illiterate.2, treat))
  X2000_infant[i] <- with(data, t.test.p(X2000_infant.2, treat))
  X2000_pop[i] <- with(data, t.test.p(X2000_pop.2, treat))
  X2000_poverty[i] <- with(data, t.test.p(X2000_poverty.2 , treat))
  p_13.1[i] <- with(data, t.test.p(p_13.2, treat))
  p_45.1[i] <- with(data, t.test.p(p_45.2, treat))
  f.nom_13.1[i] <- with(data, t.test.p(f.nom_13.2, treat))
  f.nom_45.1[i] <- with(data, t.test.p(f.nom_45.2, treat))
  g_13.1[i] <- with(data, t.test.p(g_13.2, treat))
  g_45.1[i] <- with(data, t.test.p(g_45.2, treat))
  e.nom_13.1[i] <- with(data, t.test.p(e.nom_13.2, treat))
  e.nom_45.1[i] <- with(data, t.test.p(e.nom_45.2, treat))
  Comparecimento1t.1[i] <- with(data, t.test.p(Comparecimento1t.2, treat))
  sex[i] <- with(data, t.test.p(sex.2, treat))
  age[i] <- with(data, t.test.p(age.2, treat))
}

results <-  list(
  X2000_incomepercapita,
  X2000_doctor,
  X2000_idh_education,
  X2000_idh_income,
  X2000_idh_longevity,
  X2000_illiterate,
  X2000_infant,
  X2000_pop, 
  X2000_poverty,
  p_13.1, 
  p_45.1,
  f.nom_13.1, 
  f.nom_45.1,
  g_13.1, 
  g_45.1,
  e.nom_13.1,
  e.nom_45.1,
  Comparecimento1t.1,
  sex,
  age)

pdf(paste0(dir, "paper_drafts/appendix/pvalues.pdf"),
    width=20, height=20)
par(mfrow=c(4,5), oma=c(0,3,3, 0))
for (i in 1:20){
  pv.plot(H, results[[i]], names=names[i], cexpoint=2)
}
dev.off()

#S.C.2 P-value plot, rank sum tests 

X2000_incomepercapita <- NA
X2000_doctor <- NA
X2000_idh_education <- NA
X2000_idh_income <- NA
X2000_idh_longevity <- NA
X2000_illiterate <- NA
X2000_infant <- NA
X2000_pop <- NA
X2000_poverty <- NA
p_13.1 <- NA
p_45.1 <- NA
f.nom_13.1 <- NA
f.nom_45.1 <- NA
e.nom_13.1 <- NA
e.nom_45.1 <- NA
g_13.1 <- NA
g_45.1 <- NA
Comparecimento1t.1 <- NA
age <- NA
sex <- NA

for (i in 1:length(H)){
  ## loop that first subsets our dataset 
  
  sel <-  abs(data_200315$vote_margin_share) <= H[i]
  data <- data_200315[sel,]
  control <- data[which(data$treat==0),]
  treat <- data[which(data$treat==1),]
  
  #getting the p-values
  X2000_incomepercapita[i] <- wilcox.test(treat$X2000_incomepercapita.2, 
                                          control$X2000_incomepercapita.2)$p.value
  X2000_doctor[i] <- wilcox.test(treat$X2000_doctor.2,
                                 control$X2000_doctor.2)$p.value
  X2000_idh_education[i] <- wilcox.test(treat$X2000_idh_education.2,
                                        control$X2000_idh_education.2)$p.value
  X2000_idh_income[i] <- wilcox.test(treat$X2000_idh_income.2, 
                                     control$X2000_idh_income.2)$p.value
  X2000_idh_longevity[i] <- wilcox.test(treat$X2000_idh_longevity.2, 
                                        control$X2000_idh_longevity.2)$p.value
  X2000_illiterate[i] <- wilcox.test(treat$X2000_illiterate.2, 
                                     control$X2000_illiterate.2)$p.value
  X2000_infant[i] <- wilcox.test(treat$X2000_infant.2, 
                                 control$X2000_infant.2)$p.value
  X2000_pop[i] <- wilcox.test(treat$X2000_pop.2, 
                              control$X2000_pop.2)$p.value
  X2000_poverty[i] <- wilcox.test(treat$X2000_poverty.2, 
                                  control$X2000_poverty.2)$p.value
  p_13.1[i] <- wilcox.test(treat$p_13.2, control$p_13.2)$p.value
  p_45.1[i] <- wilcox.test(treat$p_45.2, control$p_45.2)$p.value
  f.nom_13.1[i] <- wilcox.test(treat$f.nom_13.2, control$f.nom_13.2)$p.value
  f.nom_45.1[i] <- wilcox.test(treat$f.nom_45.2, control$f.nom_45.2)$p.value
  g_13.1[i] <- wilcox.test(treat$g_13.2, control$g_13.2)$p.value
  g_45.1[i] <- wilcox.test(treat$g_45.2, control$g_45.2)$p.value
  e.nom_13.1[i] <- wilcox.test(treat$e.nom_13.2, control$f.nom_13.2)$p.value
  e.nom_45.1[i] <- wilcox.test(treat$e.nom_45.2, control$e.nom_45.2)$p.value
  Comparecimento1t.1[i] <- wilcox.test(treat$Comparecimento1t.2, control$Comparecimento1t.2)$p.value
  sex[i] <- wilcox.test(treat$sex.2, control$sex.2)$p.value
  age[i] <- wilcox.test(treat$age.2, control$age.2)$p.value
}

results <-  list(
  X2000_incomepercapita,
  X2000_doctor,
  X2000_idh_education,
  X2000_idh_income,
  X2000_idh_longevity,
  X2000_illiterate,
  X2000_infant,
  X2000_pop, 
  X2000_poverty,
  p_13.1, 
  p_45.1,
  f.nom_13.1, 
  f.nom_45.1,
  g_13.1, 
  g_45.1,
  e.nom_13.1,
  e.nom_45.1,
  Comparecimento1t.1,
  sex,
  age)

pdf(paste0(dir, "paper_drafts/appendix/pvalues_wilcoxon.pdf"),
    width=20, height=20)
par(mfrow=c(4,5), oma=c(0,3,3, 0))
for (i in 1:20){
  pv.plot(H, results[[i]], names=names[i], cexpoint=2)
}
dev.off()

#S.C.3 P-value plot, F-test 

pvals <- rep(NA, length(H))
for (i in 1:length(H)){
  pvals[i] <- Fbalance(bw=H[i], data=data_200315)
}

pdf(paste0(dir, "paper_drafts/appendix/balance_ftest.pdf"), 
    height = 10, width = 10)
par(mfrow = c(1,1), oma = c(0,0, 0, 0))
pv.plot(H, pvals, name = "F-test", cexpoint = 2)
dev.off()

#S.C.4 Balance tests, local linear regression  

h <- seq(from = .01, to = .5, by = .001)

opt.bw.nsp <- optimal.bw(Y = data_200315$dtransfers.pc, 
                         X = data_200315$vote_margin_share, 
                         c = 0, reg = T)

opt.bw.mayors <- optimal.bw(Y = data_200315$dtransfers.pc.mayors, 
                            X = data_200315$vote_margin_share, 
                            c = 0, reg = T)
 
pdf(paste0(dir, "paper_drafts/appendix/balance_tests_linear.pdf"), 
    height=15, width=15)
par(mfrow=c(5,4))
for (i in 1:length(outcomes2)){
  local.linear.1(outcome=outcomes2[[i]], runvar=data_200315v4$vote_margin_share,  h, main=names[i], 
                 axes=T, opt.bw1=opt.bw.nsp, opt.bw2=opt.bw.mayors, 
                 ylab="Effect Estimates", 
                 xlab="Vote Margin (%) - Distance from Cutoff")
} 
dev.off()

#S.C.5 P-value plots, lagged outcomes (randomization inference)

H1 <- c(0.005, 0.01, 0.02, 0.03, 0.04, 0.05)

#Per capita
set.seed(06511)
est <- matrix(NA, length(H1),7)
for (i in 1:length(H1)){
  
  sel <-  abs(data_lagged$vote_margin_share) <= H1[i]
  temp <- data_lagged[sel,]     
  
  Ys <- with(temp, genouts(dtransfers.pc, treat, ate=0))
  probs <- with(temp, genprobexact(treat))
  perms <- with(temp, genperms(treat, maxiter=100000))
  ate <- with(temp, estate(dtransfers.pc, treat, prob=probs))
  distout <- with(temp, gendist(Ys, perms=perms, prob=probs))
  resultp <- dispdist(distout, ate)
  ci_ate <-  with(temp, c(invert.ci(dtransfers.pc, treat,probs,perms,0.025),
                          invert.ci(dtransfers.pc, treat,probs,perms,0.975))) 
  
  est[i,1] <- H1[i]*100
  est[i,2] <- ate
  est[i,3] <- resultp$sd
  est[i,4] <- resultp$lesser.p.value  
  est[i,5] <- resultp$two.tailed.p.value
  est[i, 6] <- ci_ate[1]
  est[i, 7] <- ci_ate[2]
}


colnames(est) <- c("Vote margin (%)", "LATE", "Standard-Deviation", 
                   "One-tailed p-value", "Two-tailed p-value", "Lower 95%", "Upper 95%")
xtable(est, digits=3)
pval <- est[, 4]
est_nsp_lagged <- est

bw <- c(0.005, 0.01, 0.02, 0.03, 0.04, 0.05)
pdf(paste0(dir, "paper_drafts/appendix/pvalues_ri_laggednsp.pdf"),
    width=8, height=8)
par(mfrow=c(1,1), oma=c(0,3,3, 0))
pv.plot(bw, est_nsp_lagged[,5], names="Lagged NSP Transfers", cexpoint=5) 
dev.off()

#Mayoral Transfers

#Per capita
est <- matrix(NA, length(H1),7)
for (i in 1:length(H1)){
  
  sel <-  abs(data_lagged$vote_margin_share) <= H1[i]
  temp <- data_lagged[sel,]     
  
  Ys <- with(temp, genouts(dtransfers.pc.mayors, treat, ate=0))
  probs <- with(temp, genprobexact(treat))
  perms <- with(temp, genperms(treat, maxiter=100000))
  ate <- with(temp, estate(dtransfers.pc.mayors, treat, prob=probs))
  distout <- with(temp, gendist(Ys, perms=perms, prob=probs))
  resultp <- dispdist(distout, ate)
  ci_ate <-  with(temp, c(invert.ci(dtransfers.pc.mayors, treat,probs,perms,0.025),
                          invert.ci(dtransfers.pc.mayors, treat,probs,perms,0.975))) 
  
  est[i,1] <- H1[i]*100
  est[i,2] <- ate
  est[i,3] <- resultp$sd
  est[i,4] <- resultp$greater.p.value  
  est[i,5] <- resultp$two.tailed.p.value
  est[i, 6] <- ci_ate[1]
  est[i, 7] <- ci_ate[2]
}

colnames(est) <- c("Vote margin (%)", "LATE", "Standard-Deviation", 
                   "One-tailed p-value", "Two-tailed p-value", "Lower 95%", "Upper 95%")
xtable(est, digits=3)
pval_mayors <- est[,3]
est_mayoral_lagged <- est

pdf(paste0(dir, "paper_drafts/appendix/pvalues_ri_laggedmayoral.pdf"),
    width=8, height=8)
par(mfrow=c(1,1), oma=c(0,3,3, 0))
pv.plot(bw, est_mayoral_lagged[,5], names="Lagged Mayoral Transfers", cexpoint=5) 
dev.off()

#S.C.6 P-value plots, lagged outcomes (normal approximation)

nsp.lagged <- NA
mayors.lagged <- NA

for (i in 1:length(H)){
  ## loop that first subsets our dataset 
  
  sel <-  abs(data_lagged$vote_margin_share) <= H[i]
  data <- data_lagged[sel,]  
  
  #getting the p-values
  nsp.lagged[i] <- with(data, t.test.p(dtransfers.pc, treat))
}  

for (i in 1:length(H)){
  ## loop that first subsets our dataset 
  
  sel <-  abs(data_lagged$vote_margin_share) <= H[i]
  data <- data_lagged[sel,]  
  
  #getting the p-values
  mayors.lagged[i] <- with(data, t.test.p(dtransfers.pc.mayors, treat))
} 

pdf(paste0(dir, "paper_drafts/appendix/pvalues_laggednsp.pdf"),
    width=8, height=8)
par(mfrow=c(1,1), oma=c(0,3,3, 0))
pv.plot(H, nsp.lagged, names="Lagged NSP Transfers", cexpoint=2)
dev.off()

pdf(paste0(dir, "paper_drafts/appendix/pvalues_laggedmayors.pdf"),
    width=8, height=8)
par(mfrow=c(1,1), oma=c(0,3,3, 0))
pv.plot(H, mayors.lagged, names="Lagged Mayoral Transfers", cexpoint=2)
dev.off()

#S.C.7 P-value plots,lagged outcomes, Wilcoxon Rank-Sum 

nsp.lagged.wilcoxon <- NA
mayors.lagged.wilcoxon <- NA

for (i in 1:length(H)){
  ## loop that first subsets our dataset 
  
  sel <-  abs(data_lagged$vote_margin_share) <= H[i]
  data <- data_lagged[sel,]
  control <- data[which(data$treat==0),]
  treat <- data[which(data$treat==1),]
  
  #getting the p-values
  nsp.lagged.wilcoxon[i] <- wilcox.test(treat$dtransfers.pc, 
                                        control$dtransfers.pc)$p.value
}

for (i in 1:length(H)){
  ## loop that first subsets our dataset 
  
  sel <-  abs(data_lagged$vote_margin_share) <= H[i]
  data <- data_lagged[sel,]
  control <- data[which(data$treat==0),]
  treat <- data[which(data$treat==1),]
  
  #getting the p-values
  mayors.lagged.wilcoxon[i] <- wilcox.test(treat$dtransfers.pc.mayors,
                                           control$dtransfers.pc.mayors)$p.value
}

pdf(paste0(dir, "paper_drafts/appendix/pvalues_wilcoxon_laggedmayoral.pdf"),
    width=8, height=8)
par(mfrow=c(1,1), oma=c(0,3,3, 0))
pv.plot(H, mayors.lagged.wilcoxon, names="Lagged Mayoral Transfers", cexpoint=2)
dev.off()

pdf(paste0(dir, "paper_drafts/appendix/pvalues_wilcoxon_laggednsp.pdf"),
    width=8, height=8)
par(mfrow=c(1,1), oma=c(0,3,3, 0))
pv.plot(H, nsp.lagged.wilcoxon, names="Lagged NSP Transfers", cexpoint=2) 
dev.off()

#S.C.8 Density of running variable (vote margin, %)

pdf(paste0(dir, "paper_drafts/appendix/density_running.pdf"),
    width=8, height=8)
mccrary <- DCdensity(data_200315$vote_margin_share, 0, ext.out=TRUE)
dev.off()

#S.C.9 RDD Plots, Full Window

pdf(paste0(dir, "paper_drafts/appendix/rdd_plot_nsp_2003_2015_cct_full.pdf"), 
    width=8, height=8)
par(mar=c(5, 8, 2.5, 0.5))
rd.plot(x = data_200315$vote_margin_share2, 
        numbinl = 100, numbinr = 100,
        y = data_200315$dtransfers.pc, y.lim = c(0, 30), x.lim= c(-100, 100),
        x.label = "President's Party Vote Margin (%)", 
        y.label = "Transfers to Non-State Providers \n (reais per capita)",
        title = "", N=FALSE)
dev.off()

pdf(paste0(dir, "paper_drafts/appendix/rdd_plot_mayoral_2003_2015_cct_full.pdf"), 
    width=8, height=8)
par(mar=c(5, 8, 2.5, 0.5))
rd.plot(x = data_200315$vote_margin_share2,
        numbinl = 150, numbinr = 150,
        y = data_200315$dtransfers.pc.mayors, y.lim = c(110, 150), x.lim=c(-100, 100),
        x.label = "President's Party Vote Margin (%)", 
        y.label = "Transfers to Mayors \n (reais per capita)", 
        title="", N=FALSE)
dev.off()

#S.C.10 External validity (Federal Transfers)

cov.05 <- as.data.frame(data_200315[abs(data_200315$vote_margin_share) <= 0.005, ])
cov.1 <- as.data.frame(data_200315[abs(data_200315$vote_margin_share) <= 0.01, ])
cov.5 <- as.data.frame(data_200315[abs(data_200315$vote_margin_share) <= 0.05, ])
cov.all <- as.data.frame(data_200315)

pdf(paste0(dir, "paper_drafts/appendix/external_all.pdf"),
    width=10, height=10)
par(mfrow=c(4,5))
m <- matrix(c(1,2,3,4,5,6,7, 8, 9, 10, 11, 12, 
              13, 14, 15, 16, 17, 18, 19, 20), 4, 5, byrow = TRUE)
layout(mat = m)
for (i in 1:20){
  par(mar = c(2,2,2,2))
  plot(density(cov.05[,86+i], bw="nrd0", adjust=.75, na.rm=T),
       lty=2, col="red", xlab = "", ylab = "Density", main = names[i])
  rug(cov.05[,86+i], ticksize = 0.03, side = 1, col="red")
  
  lines(density(cov.all[,86+i], bw="nrd0", adjust=.75, na.rm=T),
        lty=1, col="darkblue")
  rug(cov.all[,86+i], ticksize = 0.03, side = 3, col="darkblue")
}
dev.off()
#Warnings about clipping rug. It's OK.

#S.C.11 P-value plot, difference of means (state-level) 

for (i in 1:length(H)){
  ## loop that first subsets our dataset 
  
  sel <-  abs(data_200715$vote_margin_share) <= H[i]
  data <- data_200715[sel,]  
  
  #getting the p-values
  X2000_incomepercapita[i] <- with(data, t.test.p(X2000_incomepercapita.2, treat))
  X2000_doctor[i] <- with(data, t.test.p(X2000_doctor.2, treat))
  X2000_idh_education[i] <- with(data, t.test.p(X2000_idh_education.2, treat))
  X2000_idh_income[i] <- with(data, t.test.p(X2000_idh_income.2, treat))
  X2000_idh_longevity[i] <- with(data, t.test.p(X2000_idh_longevity.2, treat))
  X2000_illiterate[i] <- with(data, t.test.p(X2000_illiterate.2, treat))
  X2000_infant[i] <- with(data, t.test.p(X2000_infant.2, treat))
  X2000_pop[i] <- with(data, t.test.p(X2000_pop.2, treat))
  X2000_poverty[i] <- with(data, t.test.p(X2000_poverty.2 , treat))
  p_13.1[i] <- with(data, t.test.p(p_13.2, treat))
  p_45.1[i] <- with(data, t.test.p(p_45.2, treat))
  f.nom_13.1[i] <- with(data, t.test.p(f.nom_13.2, treat))
  f.nom_45.1[i] <- with(data, t.test.p(f.nom_45.2, treat))
  g_13.1[i] <- with(data, t.test.p(g_13.2, treat))
  g_45.1[i] <- with(data, t.test.p(g_45.2, treat))
  e.nom_13.1[i] <- with(data, t.test.p(e.nom_13.2, treat))
  e.nom_45.1[i] <- with(data, t.test.p(e.nom_45.2, treat))
  Comparecimento1t.1[i] <- with(data, t.test.p(Comparecimento1t.2, treat))
  sex[i] <- with(data, t.test.p(sex.2, treat))
  age[i] <- with(data, t.test.p(age.2, treat))
}

results <-  list(
  X2000_incomepercapita,
  X2000_doctor,
  X2000_idh_education,
  X2000_idh_income,
  X2000_idh_longevity,
  X2000_illiterate,
  X2000_infant,
  X2000_pop, 
  X2000_poverty,
  p_13.1, 
  p_45.1,
  f.nom_13.1, 
  f.nom_45.1,
  g_13.1, 
  g_45.1,
  e.nom_13.1,
  e.nom_45.1,
  Comparecimento1t.1,
  sex,
  age)

pdf(paste0(dir, "paper_drafts/appendix/pvalues_state.pdf"),
    width=20, height=20)
par(mfrow=c(4,5), oma=c(0,3,3, 0))
for (i in 1:20){
  pv.plot(H, results[[i]], names=names[i], cexpoint=2)
}
dev.off()

#S.C.12 P-value plot, rank sum tests (state-level) 

for (i in 1:length(H)){
  ## loop that first subsets our dataset 
  
  sel <-  abs(data_200715$vote_margin_share) <= H[i]
  data <- data_200715[sel,]
  control <- data[which(data$treat==0),]
  treat <- data[which(data$treat==1),]
  
  #getting the p-values
  X2000_incomepercapita[i] <- wilcox.test(treat$X2000_incomepercapita.2, 
                                          control$X2000_incomepercapita.2)$p.value
  X2000_doctor[i] <- wilcox.test(treat$X2000_doctor.2,
                                 control$X2000_doctor.2)$p.value
  X2000_idh_education[i] <- wilcox.test(treat$X2000_idh_education.2,
                                        control$X2000_idh_education.2)$p.value
  X2000_idh_income[i] <- wilcox.test(treat$X2000_idh_income.2, 
                                     control$X2000_idh_income.2)$p.value
  X2000_idh_longevity[i] <- wilcox.test(treat$X2000_idh_longevity.2, 
                                        control$X2000_idh_longevity.2)$p.value
  X2000_illiterate[i] <- wilcox.test(treat$X2000_illiterate.2, 
                                     control$X2000_illiterate.2)$p.value
  X2000_infant[i] <- wilcox.test(treat$X2000_infant.2, 
                                 control$X2000_infant.2)$p.value
  X2000_pop[i] <- wilcox.test(treat$X2000_pop.2, 
                              control$X2000_pop.2)$p.value
  X2000_poverty[i] <- wilcox.test(treat$X2000_poverty.2, 
                                  control$X2000_poverty.2)$p.value
  p_13.1[i] <- wilcox.test(treat$p_13.2, control$p_13.2)$p.value
  p_45.1[i] <- wilcox.test(treat$p_45.2, control$p_45.2)$p.value
  f.nom_13.1[i] <- wilcox.test(treat$f.nom_13.2, control$f.nom_13.2)$p.value
  f.nom_45.1[i] <- wilcox.test(treat$f.nom_45.2, control$f.nom45.2)$p.value
  g_13.1[i] <- wilcox.test(treat$g_13.2, control$g_13.2)$p.value
  g_45.1[i] <- wilcox.test(treat$g_45.2, control$g_45.2)$p.value
  e.nom_13.1[i] <- wilcox.test(treat$e.nom_13.2, control$f.nom_13.2)$p.value
  e.nom_45.1[i] <- wilcox.test(treat$e.nom_45.2, control$e.nom_45.2)$p.value
  Comparecimento1t.1[i] <- wilcox.test(treat$Comparecimento1t.2, control$Comparecimento1t.2)$p.value
  sex[i] <- wilcox.test(treat$sex.2, control$sex.2)$p.value
  age[i] <- wilcox.test(treat$age.2, control$age.2)$p.value
}

results <-  list(
  X2000_incomepercapita,
  X2000_doctor,
  X2000_idh_education,
  X2000_idh_income,
  X2000_idh_longevity,
  X2000_illiterate,
  X2000_infant,
  X2000_pop, 
  X2000_poverty,
  p_13.1, 
  p_45.1,
  f.nom_13.1, 
  f.nom_45.1,
  g_13.1, 
  g_45.1,
  e.nom_13.1,
  e.nom_45.1,
  Comparecimento1t.1,
  sex,
  age)

pdf(paste0(dir, "paper_drafts/appendix/pvalues_wilcoxon_state.pdf"),
    width=20, height=20)
par(mfrow=c(4,5), oma=c(0,3,3, 0))
for (i in 1:20){
  pv.plot(H, results[[i]], names=names[i], cexpoint=2)
}
dev.off()

#S.C.13 P-value plot, F-test (state-level) 

pvals <- rep(NA, length(H))
for (i in 1:length(H)){
  pvals[i] <- Fbalance(bw=H[i], data=data_200715)
}

#F-test
pdf(paste0(dir, "paper_drafts/appendix/balance_ftest_state.pdf"), 
    height=10, width=10)
par(mfrow=c(1,1), oma=c(0,0, 0, 0))
pv.plot(H, pvals, name="F-test", cexpoint=2)
dev.off()

#S.C.14 Balance tests, local linear regression (state-level)

opt.bw.nsp.state <- optimal.bw(Y = data_200715$dtransfers.pc, 
                         X = data_200715$vote_margin_share, 
                         c = 0, reg = T)

opt.bw.mayors.state <- optimal.bw(Y = data_200715$dtransfers.pc.mayors, 
                            X = data_200715$vote_margin_share, 
                            c = 0, reg = T)

h <- seq(from = .01, to = .5, by = .001)
pdf(paste0(dir, "paper_drafts/appendix/balance_tests_linear_state.pdf"), 
    height=15, width=15)
par(mfrow=c(5,4))
for (i in 1:length(outcomes.state)){
  local.linear.1(outcome=outcomes.state[[i]], runvar=data_200715$vote_margin_share,  h, main=names[i], 
                 axes=T, opt.bw1=opt.bw.nsp.state, opt.bw2=opt.bw.mayors.state, ylab="Effect Estimates", 
                 xlab="Vote Margin (%) - Distance from Cutoff")
} 
dev.off()

#S.C.15 Density running variable (state-level)

pdf(paste0(dir, "paper_drafts/appendix/state_density_running.pdf"),
    width=8, height=8)
mccrary <- DCdensity(data_200715$vote_margin_share, 0, ext.out=TRUE)
dev.off()

#S.C.16 State NSP and Mayoral Transfers: RDD Estimates, 2007–2015 

pdf(paste0(dir, "paper_drafts/appendix/nsp_rdd_disbursed_perc_2007_2015_state.pdf"),
    width=8, height=8)
RD_plot(running=data_200715$vote_margin_share2, treat=data_200715$treat, 
        outcome=data_200715$dtransfers.pc, clusterT=T, 
        cluster.var=data_200715$ibge2,
        cutoff=0, min_running=0.2, max_running=5,
        opt_bw=F, nr_windows=20, main_est="dom", ci="90%",
        add_est=F, nr_obs=T, 
        nr_obs_lab=c(10, 25, 50, 100, 150),
        label_x="Vote margin (%)\n Absolute Distance from Cutoff",
        label_y="Average Difference in NSP Transfers \n Between Aligned and Unaligned",
        plot_label="",
        legend=T)
dev.off()

pdf(paste0(dir, "paper_drafts/appendix/mayors_rdd_disbursed_perc_2007_2015_state.pdf"),
    width=8, height=8)
RD_plot(running=data_200715$vote_margin_share2, treat=data_200715$treat, 
        outcome=data_200715$dtransfers.pc.mayors, clusterT=T, 
        cluster.var=data_200715$ibge2,
        cutoff=0, min_running=0.2, max_running=5,
        opt_bw=F, nr_windows=20, main_est="dom", ci="90%",
        add_est=F, nr_obs=T, 
        nr_obs_lab=c(10, 25, 50, 100, 150),
        label_x="Vote margin (%)\n Absolute Distance from Cutoff",
        label_y="Average Difference in NSP Transfers \n Between Aligned and Unaligned",
        plot_label="",
        legend=T)
dev.off()


#S.C.17 External validity (State Transfers)

cov.05 <- as.data.frame(data_200715[abs(data_200715$vote_margin_share) <= 0.005, ])
cov.1 <- as.data.frame(data_200715[abs(data_200715$vote_margin_share) <= 0.01, ])
cov.5 <- as.data.frame(data_200715[abs(data_200715$vote_margin_share) <= 0.05, ])
cov.all <- as.data.frame(data_200715)

pdf(paste0(dir, "paper_drafts/appendix/external_all_state.pdf"),
    width=10, height=10)
par(mfrow=c(4,5))
m <- matrix(c(1,2,3,4,5,6,7, 8, 9, 10, 11, 12, 
              13, 14, 15, 16, 17, 18, 19, 20), 4, 5, byrow = TRUE)
layout(mat = m)
for (i in 1:20){
  par(mar = c(2,2,2,2))
  plot(density(cov.05[,82+i], bw="nrd0", adjust=.75, na.rm=T),
       lty=2, col="red", xlab = "", ylab = "Density", main = names[i])
  rug(cov.05[,82+i], ticksize = 0.03, side = 1, col="red")
  
  lines(density(cov.all[,82+i], bw="nrd0", adjust=.75, na.rm=T),
        lty=1, col="darkblue")
  rug(cov.all[,82+i], ticksize = 0.03, side = 3, col="darkblue")
}
dev.off()
#Warnings about clipping rug. It's OK.

#S.C.18 Legislative Amendments, Difference in Means, NSP and Mayoral Transfers,
# 2003–2011 (in July 2012 reais)

#Cleaning
emendasf <- emendasf[!is.na(emendasf$IBGECod),]

pdf(paste0(dir, "paper_drafts/appendix/nsp_rdd_disbursed_perc_2002_2011.pdf"),
    width=8, height=8)
RD_plot(running = emendasf$vote_margin_share2, treat = emendasf$treat, 
        outcome = emendasf$dliquidadosum0.pc, 
        cutoff = 0, min_running = 0.2, max_running = 5, clusterT = T,
        cluster.var = emendasf$IBGECod,
        opt_bw = F, nr_windows = 20, main_est="dom", ci="90%",
        add_est = F, nr_obs = T, 
        nr_obs_lab = c(30, 250, 500, 700, 900),
        label_x = "Vote margin (%)\n Absolute Distance from Cutoff",
        label_y = "Average Difference in NSP Transfers \n Between Aligned and Unaligned",
        plot_label = "",
        legend = T)
dev.off()

pdf(paste0(dir, "paper_drafts/appendix/mayor_rdd_disbursed_pc_2002_2011.pdf"),
    width=8, height=8)
RD_plot(running = emendasf$vote_margin_share2, treat = emendasf$treat, 
        outcome = emendasf$dliquidadosum0.pc.mayors, 
        cutoff = 0, min_running = 0.2, max_running = 5,
        clusterT = T,
        cluster.var = emendasf$IBGECod,
        opt_bw = F, nr_windows = 20, main_est = "dom", ci = "90%",
        add_est = F, nr_obs = T, 
        nr_obs_lab = c(30, 250, 500, 700, 900),
        label_x = "Vote margin (%)\n Absolute Distance from Cutoff",
        label_y = "Average Difference in Mayoral Transfers \n Between Aligned and Unaligned",
        plot_label = "",
        legend = T)
dev.off()

#S.D.1 Carryover Effects

task1 <- amce(d_choice ~ nsp_cityff + federalff + partisanshipff,
                 data=conjointff[conjointff$task==1,],
                 cluster=TRUE, respondent.id="respID")

task2 <- amce(d_choice ~ nsp_cityff + federalff + partisanshipff,
             data=conjointff[conjointff$task==2,],
             cluster=TRUE, respondent.id="respID")

task3 <- amce(d_choice ~ nsp_cityff + federalff + partisanshipff,
                data=conjointff[conjointff$task==3,],
                cluster=TRUE, respondent.id="respID")

task4 <- amce(d_choice ~ nsp_cityff + federalff + partisanshipff,
              data=conjointff[conjointff$task==4,],
              cluster=TRUE, respondent.id="respID")

pdf(paste0(dir, "paper_drafts/appendix/carryover1.pdf"), width=7, height=7)
par(mar=c(4.1,3.1,4.1,1.1))
plot_cjoint(results = task1)
dev.off()

pdf(paste0(dir, "paper_drafts/appendix/carryover2.pdf"), width=7, height=7)
par(mar=c(4.1,3.1,4.1,1.1))
plot_cjoint(results = task2)
dev.off()

pdf(paste0(dir, "paper_drafts/appendix/carryover3.pdf"), width=7, height=7)
par(mar=c(4.1,3.1,4.1,1.1))
plot_cjoint(results = task3)
dev.off()

pdf(paste0(dir, "paper_drafts/appendix/carryover4.pdf"), width=7, height=7)
par(mar=c(4.1,3.1,4.1,1.1))
plot_cjoint(results = task4)
dev.off()

#S.D.2 Profile Order Effects

profile1 <- amce(d_choice ~ nsp_cityff + federalff + partisanshipff,
              data=conjointff[conjointff$panel=="A",],
              cluster=TRUE, respondent.id="respID")

profile2 <- amce(d_choice ~ nsp_cityff + federalff + partisanshipff,
              data=conjointff[conjointff$panel=="B",],
              cluster=TRUE, respondent.id="respID")

pdf(paste0(dir, "paper_drafts/appendix/panel1.pdf"), width=7, height=7)
par(mar=c(4.1,3.1,4.1,1.1))
plot_cjoint(results = profile1)
dev.off()

pdf(paste0(dir, "paper_drafts/appendix/panel2.pdf"), width=7, height=7)
par(mar=c(4.1,3.1,4.1,1.1))
plot_cjoint(results = profile2)
dev.off()

#S.D.3 Federal Government Outcome Results

#Rate Outcome 
results.rate.fed.psdb <- amce(ratefedf ~ nsp_cityff + federalff,
                              data=conjointff.psdb,
                              cluster=TRUE, respondent.id="respID")
results.rate.fed.pt <- amce(ratefedf ~ nsp_cityff + federalff,
                            data=conjointff.pt,
                            cluster=TRUE, respondent.id="respID")


pdf(paste0(dir, "paper_drafts/appendix/federal_rate_unaligned.pdf"), width=7, height=7)
par(mar=c(4.1,3.1,4.1,1.1))
plot_cjoint(results = results.rate.fed.psdb)
dev.off()

pdf(paste0(dir, "paper_drafts/appendix/federal_rate_aligned.pdf"), width=7, height=7)
par(mar=c(4.1,3.1,4.1,1.1))
plot_cjoint(results = results.rate.fed.pt)
dev.off()  

#S.D.4 Effects of Welfare and Joint Provision on Probability of Mayor being Preferred
#and Change of Mayoral Ratings (weigthed survey sample)

#Recoding for analysis
conjointff$nsp_cityf <- as.factor(conjointff$nsp_cityf)
conjointff$federalf <- as.factor(conjointff$federalf)

#Analysis by aligned and unaligned
conjointff.psdb <- conjointff %>% filter(partisanship == "O prefeito pertence ao PSDB (número 45).")
conjointff.pt <- conjointff %>% filter(partisanship == "O prefeito pertence ao  PT (número 13).")
conjointff.psdb$nsp_cityf <- relevel(conjointff.psdb$nsp_cityf, ref="A PREFEITURA é responsável pelo abrigo para moradores de rua.")
conjointff.pt$nsp_cityf <- relevel(conjointff.pt$nsp_cityf, ref="A PREFEITURA é responsável pelo abrigo para moradores de rua.")

#Choice
results.choice.psdbw <- amce(d_choice ~ nsp_cityf + federalf,
                            data=conjointff.psdb,
                            cluster=TRUE, respondent.id="respID", weights="weights")
results.choice.ptw <- amce(d_choice ~ nsp_cityf + federalf,
                          data=conjointff.pt,
                          cluster=TRUE, respondent.id="respID", weights="weights")

#Rate
results.rate.psdbw <- amce(ratef ~ nsp_cityf + federalf,
                          data=conjointff.psdb,
                          cluster=TRUE, respondent.id="respID", weights="weights")
results.rate.ptw <- amce(ratef ~ nsp_cityf + federalf,
                        data=conjointff.pt,
                        cluster=TRUE, respondent.id="respID", weights="weights")

pdf(paste0(dir, "paper_drafts/appendix/weigthed_choicepsdb.pdf"), width=7, height=7)
par(mar=c(4.1,3.1,4.1,1.1))
plot_cjoint(results = results.choice.psdbw)
dev.off()

pdf(paste0(dir, "paper_drafts/appendix/weigthed_choicept.pdf"), width=7, height=7)
par(mar=c(4.1,3.1,4.1,1.1))
plot_cjoint(results = results.choice.ptw)
dev.off()

pdf(paste0(dir, "paper_drafts/appendix/weigthed_ratepsdb.pdf"), width=7, height=7)
par(mar=c(4.1,3.1,4.1,1.1))
plot_cjoint(results = results.rate.psdbw, x_axis = c(-2, 2))
dev.off()

pdf(paste0(dir, "paper_drafts/appendix/weigthed_ratept.pdf"), width=7, height=7)
par(mar=c(4.1,3.1,4.1,1.1))
plot_cjoint(results = results.rate.ptw, x_axis = c(-2, 2))
dev.off()

length(!is.na(conjointff.psdb$weights))
length(!is.na(conjointff.pt$weights))

#S.D.5 Effects of Welfare and Joint Provision on Probability of Mayor being Preferred
#and Change of Mayoral Ratings (lower income)

#recoding
conjointffl.psdb <- conjointff.psdb[conjointff.psdb$MWf==1,]
conjointffl.pt <- conjointff.pt[conjointff.pt$MWf==1,]

#Choice
results.choice.psdbl <- amce(d_choice ~ nsp_cityf + federalf,
                             data=conjointffl.psdb,
                             cluster=TRUE, respondent.id="respID")
results.choice.ptl <- amce(d_choice ~ nsp_cityf + federalf,
                           data=conjointffl.pt,
                           cluster=TRUE, respondent.id="respID")

#Rate
results.rate.psdbl <- amce(ratef ~ nsp_cityf + federalf,
                           data=conjointffl.psdb,
                           cluster=TRUE, respondent.id="respID")
results.rate.ptl <- amce(ratef ~ nsp_cityf + federalf,
                         data=conjointffl.pt,
                         cluster=TRUE, respondent.id="respID")


pdf(paste0(dir, "paper_drafts/appendix/lower_choicept.pdf"),
    width=7, height=7)
par(mar=c(4.1,3.1,4.1,1.1))
plot_cjoint(results = results.choice.ptl)
dev.off()

pdf(paste0(dir, "paper_drafts/appendix/lower_choicepsdb.pdf"),
    width=7, height=7)
par(mar=c(4.1,3.1,4.1,1.1))
plot_cjoint(results = results.choice.psdbl)
dev.off()

pdf(paste0(dir, "paper_drafts/appendix/lower_ratept.pdf"),
    width=7, height=7)
par(mar=c(4.1,3.1,4.1,1.1))
plot_cjoint(results = results.rate.ptl, x_axis = c(-2, 2))
dev.off()

pdf(paste0(dir, "paper_drafts/appendix/lower_ratepsdb.pdf"),
    width=7, height=7)
par(mar=c(4.1,3.1,4.1,1.1))
plot_cjoint(results = results.rate.psdbl, x_axis = c(-2, 2))
dev.off()

conjointffl <- bind_rows(conjointffl.pt, conjointffl.psdb)

length(unique(conjointffl$respID)) 
length(!is.na(conjointffl.psdb$d_choice))
length(!is.na(conjointffl.pt$d_choice))

#S.D.6 Effects of Welfare and Joint Provision on Probability of Mayor being Preferred
#and Changeof Mayoral Ratings (co-partisans)

#Analysis
conjointff$copartisan <- ifelse(((conjointff$partisanship=="O prefeito pertence ao PSDB (número 45)." &
                                    conjointff$party=="PSDB (Partido da Social Demoracia Brasileira)") |
                                   (conjointff$partisanship=="O prefeito pertence ao  PT (número 13)." &
                                      conjointff$party=="PT (Partido dos Trabalhadores)")), 1, 0)


conjointff$copartisanvote <- ifelse(((conjointff$partisanship=="O prefeito pertence ao PSDB (número 45)." &
                                        conjointff$vote2014=="Aécio Neves (PSDB, 45)") |
                                       (conjointff$partisanship=="O prefeito pertence ao  PT (número 13)." &
                                          conjointff$vote2014=="Dilma Rousseff (PT, 13)")), 1, 0)

conjointffc <- conjointff[conjointff$copartisanvote==1,]
conjointffc.psdb <- conjointffc[conjointffc$partisanship=="O prefeito pertence ao PSDB (número 45).",]
conjointffc.pt <- conjointffc[conjointffc$partisanship=="O prefeito pertence ao  PT (número 13).",]

#Choice
results.choice.psdbc <- amce(d_choice ~ nsp_cityf + federalf,
                             data=conjointffc.psdb,
                             cluster=TRUE, respondent.id="respID")
results.choice.ptc <- amce(d_choice ~ nsp_cityf + federalf,
                           data=conjointffc.pt,
                           cluster=TRUE, respondent.id="respID")

#Rate
results.rate.psdbc <- amce(ratef ~ nsp_cityf + federalf,
                           data=conjointffc.psdb,
                           cluster=TRUE, respondent.id="respID")
results.rate.ptc <- amce(ratef ~ nsp_cityf + federalf,
                         data=conjointffc.pt,
                         cluster=TRUE, respondent.id="respID")


pdf(paste0(dir, "paper_drafts/appendix/co_choicept.pdf"),
    width=7, height=7)
par(mar=c(4.1,3.1,4.1,1.1))
plot_cjoint(results = results.choice.ptc)
dev.off()

pdf(paste0(dir, "paper_drafts/appendix/co_choicepsdb.pdf"),
    width=7, height=7)
par(mar=c(4.1,3.1,4.1,1.1))
plot_cjoint(results = results.choice.psdbc)
dev.off()

pdf(paste0(dir, "paper_drafts/appendix/co_ratept.pdf"),
    width=7, height=7)
par(mar=c(4.1,3.1,4.1,1.1))
plot_cjoint(results = results.rate.ptc, x_axis = c(-2, 2))
dev.off()

pdf(paste0(dir, "paper_drafts/appendix/co_ratepsdb.pdf"),
    width=7, height=7)
par(mar=c(4.1,3.1,4.1,1.1))
plot_cjoint(results = results.rate.psdbc, x_axis = c(-2, 2))
dev.off()

length(unique(conjointffc$respID))
length(!is.na(conjointffc.psdb$weights))
length(!is.na(conjointffc.pt$weights))

########## Data cited in the paper and in the Appendix

#Section 2.2

#Distribution of Transfers
type_total <- sic %>% group_by(TX_ESFERA_ADM_PROPONENTE) %>%
			summarise(sum_type = sum(d_empenhado)) 
(type_total$sum_type/sum(type_total$sum_type))*100			
			
#Programs and Type of Provider
sic <- sic %>% mutate(nsp = ifelse(TX_ESFERA_ADM_PROPONENTE == "PRIVADA", 1, ifelse(TX_ESFERA_ADM_PROPONENTE == "MUNICIPAL", 0, NA)))

programs <- table(sic$CD_PROGRAMA, sic$nsp)
bothPrograms <- programs[(programs[,1]>0 & programs[,2]>0),]
muns_bothp <- unique(rownames(bothPrograms))

acao <- table(sic$CD_ACAO_PPA, sic$nsp)
bothAcao <- acao[(acao[,1]>0 & acao[,2]>0),]
muns_botha <- unique(rownames(bothAcao))

sic$acao_both[sic$CD_ACAO_PPA %in% muns_botha] <- 1
sic$acao_both[!(sic$CD_ACAO_PPA %in% muns_botha)] <- 0

sic$program_both[sic$CD_PROGRAMA %in% muns_bothp] <- 1
sic$program_both[!(sic$CD_PROGRAMA %in% muns_bothp)] <- 0

prop.table(table(sic$program_both))
prop.table(table(sic$acao_both))

#Section 3.2		

summary(data_200315$bf_pc)
summary(data_200315$iptu_pc)
summary(data_200315$as_prev_pc)
	
#Revoking or boosting transfers?			
#NSP
r2000 <- all_close[[1]] 
r2004 <- all_close[[2]]
r2008 <- all_close[[3]]

r1 <- apply(r2000, 2, mean, na.rm=T)
r2 <- apply(r2004, 2, mean, na.rm=T)
r3 <- apply(r2008, 2, mean, na.rm=T)

#Combined effect
#previous year - later year

as.numeric(r1[20]) + as.numeric(r2[21]) + as.numeric(r3[22])
#It should be positive for NSPs because it's (t-1)_npt - t_pt

#Mayors

r2000 <- all_closemayors[[1]]
r2004 <- all_closemayors[[2]]
r2008 <- all_closemayors[[3]]

r1 <- apply(r2000, 2, mean, na.rm=T)
r2 <- apply(r2004, 2, mean, na.rm=T)
r3 <- apply(r2008, 2, mean, na.rm=T)

#Combined effect
#previous year (transfers) - later year (transfers)

as.numeric(r1[20]) + as.numeric(r2[21]) + as.numeric(r3[22])
#It should be negative for mayors because it's (t-1)_npt - t_pt

#Section 4.1

#Are NSP provide with federal funds stat. different from state provider with federal funds?
state_fed <- 0.027
se_state_fed <- 0.02
nsp_fed <- 0.116
se_nsp_fed <- 0.021
se_diff <- sqrt((se_state_fed^2) + (se_nsp_fed^2))
(nsp_fed - state_fed)/se_diff
			  
#Section 4.3

prop.table(table(is.na(orgs_data$party_director)))

#Boards
all_boardscnes <- all_boardscnes %>% mutate(party = ifelse(!is.na(NOME.DO.PARTIDO), 1, 0))

boards <- aggregate(party ~ cnpj_clean, sum, data=all_boardscnes)
nrow(boards[boards$party > 0,])/nrow(boards) 

ind <- aggregate(ind ~ cnpj_clean, sum, data=all_boardscnes)
all <- cbind(boards, ind)
mean(all[,2]/all[,4])
summary(all[,2]/all[,4])

prop.table(table(all_boardscnes$NOME.DO.PARTIDO))*100

#Party dominance
pr <- table(all_boardscnes$cnpj_clean, all_boardscnes$NOME.DO.PARTIDO)
tab <- cbind(pr, Total = rowSums(pr))
summary(tab)
r <- table(tab[,ncol(tab)])
14/sum(r[-1])

#Party switchers
pp <- table(all_boardscnes$NUMERO.DA.INSCRICAO, all_boardscnes$NOME.DO.PARTIDO)
tab <- cbind(pp, Total = rowSums(pp))
summary(tab)
r <- table(tab[,ncol(tab)])
240/sum(r[-1])

#Section 5 

table(cepim_orgs$receipt_closerace)

########### Additional analyses not shown in appendix

# Additional robustness
# Removing 0, -100 and 100 in vote margin

data_200315_r <- data_200315 %>% filter(vote_margin_share2 !=0, vote_margin_share2 != 100, vote_margin_share != -100)

#Optimal bandwidths
nsp_cctbw <- rdbwselect(y = data_200315_r$dtransfers.pc,  
							x= data_200315_r$vote_margin_share,
                        all = TRUE, kernel = "uniform", p=  1, q = 2)
                        
mayors_cctbw <- rdbwselect(y = data_200315_r$dtransfers.pc.mayors, 
								x = data_200315_r$vote_margin_share,
                           all = TRUE, kernel = "tri", p = 1, q = 2)

opt.bw.nsp <- optimal.bw(Y = data_200315_r$dtransfers.pc, 
                         X = data_200315_r$vote_margin_share, 
                         c = 0, reg = T)

opt.bw.mayors <- optimal.bw(Y = data_200315_r$dtransfers.pc.mayors, 
                            X = data_200315_r$vote_margin_share, 
                            c = 0, reg = T)

#NSP
rdd.table(runvar = data_200315_r$vote_margin_share, 
		  treat = data_200315_r$treat,
          outcome.1 = data_200315_r$dtransfers.pc, clusterT = T, 
          cluster.var = data_200315_r$ibge2,
          non.zero = data_200315_r$nonzero,
           p = 1, q = 2, kernel = "uniform",
          opt.pc = opt.bw.nsp)

#Mayors
rdd.table(runvar = data_200315_r$vote_margin_share, 
		  treat= data_200315_r$treat,
          outcome.1 = data_200315_r$dtransfers.pc.mayors, 
          clusterT = T, cluster.var=data_200315_r$ibge2,
          non.zero = data_200315_r$nonzero.mayors, 
          p = 1, q = 2, kernel = "tri",
          opt.pc = opt.bw.mayors)

# Controls 
#Election dummies for controls
data_200315$e2000 <- ifelse(data_200315$ANO_ELEICAO == 2000, 1, 0)
data_200315$e2004 <- ifelse(data_200315$ANO_ELEICAO == 2004, 1, 0)
data_200315$e2008 <- ifelse(data_200315$ANO_ELEICAO == 2008, 1, 0)
data_200315$e2012 <- ifelse(data_200315$ANO_ELEICAO == 2012, 1, 0)
data_200315$ANO_ELEICAO2 <- as.numeric(as.factor(data_200315$ANO_ELEICAO))

summary(rdrobust(y = data_200315$dtransfers.pc, x = data_200315$vote_margin_share, 
                 cluster=data_200315$ibge2, covs = data_200315$age.2,
                 p = 1, q = 2, kernel = "uniform"))

summary(rdrobust(y = data_200315$dtransfers.pc, x = data_200315$vote_margin_share, 
                 cluster=data_200315$ibge2, covs = data_200315$f.nom_45.2,
                 p = 1, q = 2, kernel = "uniform"))
                 
summary(rdrobust(y = data_200315$dtransfers.pc, x = data_200315$vote_margin_share, 
                 cluster=data_200315$ibge2, covs = data_200315$X2000_doctor.2,
                 p = 1, q = 2, kernel = "uniform")) 
                 
summary(rdrobust(y = data_200315$dtransfers.pc, x = data_200315$vote_margin_share, 
                 cluster=data_200315$ibge2, covs = data_200315$X2000_idh_longevity.2,
                 p = 1, q = 2, kernel = "uniform"))                
                 
summary(rdrobust(y = data_200315$dtransfers.pc, x = data_200315$vote_margin_share, 
                 cluster=data_200315$ibge2, covs = data_200315$e.nom_45.2,
                 p = 1, q = 2, kernel = "uniform"))    
                 
summary(rdrobust(y = data_200315$dtransfers.pc, x = data_200315$vote_margin_share, 
                 cluster=data_200315$ibge2, covs = data_200315$Comparecimento1t.2,
                 p = 1, q = 2, kernel = "uniform"))                                            

summary(rdrobust(y = data_200315$dtransfers.pc, x = data_200315$vote_margin_share, 
                 cluster=data_200315$ibge2, covs = data_200315$ANO_ELEICAO2,
                 p = 1, q = 2, kernel = "uniform"))

summary(rdrobust(y = data_200315$dtransfers.pc, x = data_200315$vote_margin_share, 
                 cluster=data_200315$ibge2, covs = data_200315$e2000,
                 p = 1, q = 2, kernel = "uniform"))

summary(rdrobust(y = data_200315$dtransfers.pc, x = data_200315$vote_margin_share, 
                 cluster=data_200315$ibge2, covs = data_200315$e2004,
                 p = 1, q = 2, kernel = "uniform"))

summary(rdrobust(y = data_200315$dtransfers.pc, x = data_200315$vote_margin_share, 
                 cluster=data_200315$ibge2, covs = data_200315$e2008,
                 p = 1, q = 2, kernel = "uniform"))

summary(rdrobust(y = data_200315$dtransfers.pc, x = data_200315$vote_margin_share, 
                 cluster=data_200315$ibge2, covs = data_200315$e2012,
                 p = 1, q = 2, kernel = "uniform"))
                 
                
summary(rdrobust(y = data_200315$dtransfers.pc.mayors, x = data_200315$vote_margin_share, 
                 cluster=data_200315$ibge2, covs = data_200315$age.2,
                 p = 1, q = 2, kernel = "tri"))

summary(rdrobust(y = data_200315$dtransfers.pc.mayors, x = data_200315$vote_margin_share, 
                 cluster=data_200315$ibge2, covs =data_200315$f.nom_45.2,
                 p = 1, q = 2, kernel = "tri"))

summary(rdrobust(y = data_200315$dtransfers.pc.mayors, x = data_200315$vote_margin_share, 
                 cluster=data_200315$ibge2, covs =data_200315$X2000_doctor.2,
                 p = 1, q = 2, kernel = "tri"))
                 
summary(rdrobust(y = data_200315$dtransfers.pc.mayors, x = data_200315$vote_margin_share, 
                 cluster=data_200315$ibge2, covs =data_200315$e.nom_45.2,
                 p = 1, q = 2, kernel = "tri"))
                 
 summary(rdrobust(y = data_200315$dtransfers.pc.mayors, x = data_200315$vote_margin_share, 
                 cluster=data_200315$ibge2, covs =data_200315$X2000_idh_longevity.2,
                 p = 1, q = 2, kernel = "tri"))
                 
summary(rdrobust(y = data_200315$dtransfers.pc.mayors, x = data_200315$vote_margin_share, 
                 cluster=data_200315$ibge2, covs =data_200315$Comparecimento1t.2,
                 p = 1, q = 2, kernel = "tri"))    

summary(rdrobust(y = data_200315$dtransfers.pc.mayors, x = data_200315$vote_margin_share, 
                 cluster=data_200315$ibge2, covs = data_200315$ANO_ELEICAO2,
                 p = 1, q = 2, kernel = "tri"))

summary(rdrobust(y = data_200315$dtransfers.pc.mayors, x = data_200315$vote_margin_share, 
                 cluster=data_200315$ibge2, covs = data_200315$e2000, 
                 p = 1, q = 2, kernel = "tri"))

summary(rdrobust(y = data_200315$dtransfers.pc.mayors, x = data_200315$vote_margin_share, 
                 cluster=data_200315$ibge2, covs = data_200315$e2004,
                 p = 1, q = 2, kernel = "tri"))

summary(rdrobust(y = data_200315$dtransfers.pc.mayors, x = data_200315$vote_margin_share, 
                 cluster=data_200315$ibge2, covs = data_200315$e2008,
                 p = 1, q = 2, kernel = "tri"))

summary(rdrobust(y = data_200315$dtransfers.pc.mayors, x = data_200315$vote_margin_share, 
                 cluster=data_200315$ibge2, covs = data_200315$e2012,
                 p = 1, q = 2, kernel = "tri"))          

#Using different kernels in main specifications

#Optimal bandwidths
nsp_cctbw <- rdbwselect(y = data_200315$dtransfers.pc,  
							x= data_200315$vote_margin_share,
                        all = TRUE, kernel = "tri", p=  1, q = 2)
                        
mayors_cctbw <- rdbwselect(y = data_200315$dtransfers.pc.mayors, 
								x = data_200315$vote_margin_share,
                           all = TRUE, kernel = "uniform", p = 1, q = 2)

opt.bw.nsp <- optimal.bw(Y = data_200315$dtransfers.pc, 
                         X = data_200315$vote_margin_share, 
                         c = 0, reg = T)

opt.bw.mayors <- optimal.bw(Y = data_200315$dtransfers.pc.mayors, 
                            X = data_200315$vote_margin_share, 
                            c = 0, reg = T)

#NSP
rdd.table(runvar = data_200315$vote_margin_share, treat = data_200315$treat,
          outcome.1 = data_200315$dtransfers.pc, clusterT = T, 
          cluster.var = data_200315$ibge2,
          non.zero = data_200315$nonzero, p = 1, q = 2, kernel = "tri",
          opt.pc = opt.bw.nsp)

#Mayors
rdd.table(runvar = data_200315$vote_margin_share, treat= data_200315$treat,
          outcome.1 = data_200315$dtransfers.pc.mayors, 
          clusterT = T, cluster.var=data_200315$ibge2,
          non.zero = data_200315$nonzero.mayors, p = 1, q = 2, kernel = "uniform",
          opt.pc = opt.bw.mayors)

## Excluding zeros federal transfers

data_200315v2 <- data_200315 %>% filter(dtransfers.pc > 0)
data_200315v3 <- data_200315 %>% filter(dtransfers.pc.mayors > 0)

#Optimal bandwidths
nsp_cctbw <- rdbwselect(y = data_200315v2$dtransfers.pc,  
                        x= data_200315v2$vote_margin_share,
                        all = TRUE, kernel = "uniform", p=  1, q = 2)

mayors_cctbw <- rdbwselect(y = data_200315v3$dtransfers.pc.mayors, 
                           x = data_200315v3$vote_margin_share,
                           all = TRUE, kernel = "tri", p = 1, q = 2)

opt.bw.nsp <- optimal.bw(Y = data_200315v2$dtransfers.pc, 
                         X = data_200315v2$vote_margin_share, 
                         c = 0, reg = T)

opt.bw.mayors <- optimal.bw(Y = data_200315v3$dtransfers.pc.mayors, 
                            X = data_200315v3$vote_margin_share, 
                            c = 0, reg = T)

#NSP
rdd.table(runvar = data_200315v2$vote_margin_share, treat = data_200315v2$treat,
          outcome.1 = data_200315v2$dtransfers.pc, clusterT = T, 
          cluster.var = data_200315v2$ibge2,
          non.zero = data_200315v2$nonzero, p = 1, q = 2, kernel = "uniform",
          opt.pc = opt.bw.nsp)

#Mayors
rdd.table(runvar = data_200315v3$vote_margin_share, treat= data_200315v3$treat,
          outcome.1 = data_200315v3$dtransfers.pc.mayors, 
          clusterT = T, cluster.var=data_200315v3$ibge2,
          non.zero = data_200315v3$nonzero.mayors, p = 1, q = 2, kernel = "tri",
          opt.pc = opt.bw.mayors)

#Excluding zeros state transfers

data_200715v2 <- data_200715 %>% filter(dtransfers.pc > 0)
data_200715v3 <- data_200715 %>% filter(dtransfers.pc.mayors > 0)

#Optimal bandwidths
nsp_cctbw <- rdbwselect(y = data_200715v2$dtransfers.pc,  
                        x= data_200715v2$vote_margin_share,
                        all = TRUE, kernel = "uniform", p=  1, q = 2)

mayors_cctbw <- rdbwselect(y = data_200715v3$dtransfers.pc.mayors, 
                           x = data_200715v3$vote_margin_share,
                           all = TRUE, kernel = "tri", p = 1, q = 2)

opt.bw.nsp <- optimal.bw(Y = data_200715v2$dtransfers.pc, 
                         X = data_200715v2$vote_margin_share, 
                         c = 0, reg = T)

opt.bw.mayors <- optimal.bw(Y = data_200715v3$dtransfers.pc.mayors, 
                            X = data_200715v3$vote_margin_share, 
                            c = 0, reg = T)

#NSP
rdd.table(runvar = data_200715v2$vote_margin_share, treat = data_200715v2$treat,
          outcome.1 = data_200715v2$dtransfers.pc, clusterT = T, 
          cluster.var = data_200715v2$ibge2,
          non.zero = data_200715v2$nonzero, p = 1, q = 2, kernel = "uniform",
          opt.pc = opt.bw.nsp)

#Mayors
rdd.table(runvar = data_200715v3$vote_margin_share, treat= data_200715v3$treat,
          outcome.1 = data_200715v3$dtransfers.pc.mayors, 
          clusterT = T, cluster.var=data_200715v3$ibge2,
          non.zero = data_200715v3$nonzero.mayors, p = 1, q = 2, kernel = "tri",
          opt.pc = opt.bw.mayors)

#Attribute order effects
conjointff <- conjointff %>% mutate(row1f = as.factor(row1), 
                                    row2f = as.factor(row2), 
                                    row3f = as.factor(row3))

rows <- amce(d_choice ~ nsp_cityff + federalff + partisanshipff 
              + row1f + row2f + row3f
              +  nsp_cityff*row1f + nsp_cityff*row2f +   nsp_cityff*row3f
              + federalff*row1f + federalff*row2f +   federalff*row3f
              + partisanshipff*row1f + partisanshipff*row2f +  partisanshipff*row3f,
              data=conjointff,
              cluster=TRUE, respondent.id="respID")

plot_cjoint(results = rows)
